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1 Significance of the Problem 


It has been long recognized that airborne sound transmission in the human auditory system 
(through the ear canal to the middle ear and to the inner ear) is not the only physical mech¬ 
anism responsible for hearing. The other mechanisms (generally termed “bone conduction”, 
although not only the bone may be involved) include vibrations of the ear canal walls, direct 
mechanical excitation of the ossicular system and of the cochlea through bone and soft tissue 
vibrations, excitation of the middle ear by waves passing through the nose cavity, and similar 
phenomena. While those mechanisms may be beneficial in the context of hearing aids, in the 
noisy environment they constitute a significant risk, which may be difficult to mitigate by means 
of protective devices, typically designed to reduce the airborne energy transfer. Numerical sim¬ 
ulation of energy transfer processes may be, therefore, of significant value in understanding the 
risks and preventing them by an appropriate design of protective devices, especially considering 
the fact that many relevant physical quantities are difficult and sometimes impossible to mea¬ 
sure experimentally (particularly in vivo), and physiological measurements (such as hearing 
thresholds) may be biased by individual subjective perception differences. 

2 Project Objectives 

The objectives of this project were: 

■ to develop a formulation capable of describing energy flow in a human head subject to an 
incident acoustic wave propagating in the surrounding air, and construct a corresponding 
accurate, error controlled, and efficient numerical simulation tool, 

and, using the developed simulation tool, 

■ to investigate and assess quantitatively various mechanisms of energy transfer to the 
inner ear due to various sound-wave conduction paths. 


3 Overview of the approach 

The present project consisted of three main parts: 

■ development of a fast surface integral equation solver with FFT-based matrix compres¬ 
sion, 

■ development of an anatomically faithful geometry model 

- containing all essential anatomical elements needed to study energy transfer to the 
cochlea region, and 

- fulfilling all the requirements ensuring its usability in numerical simulations. 

■ analysis employing the developed geometry and the new solver which had as a goal 
providing information on energy flow to the inner ear area. 
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A very unique and non-trivial algorithmic challenge of constructing an accurate and efficient 
solution procedure for the problem of a human body/human head surrounded by air and 
exposed to the incoming acoustic wave is the high contrast nature of the problem. Due to the 
large (of the order of 10 3 ) density ratio of human versus air tissues there is a large impedance 
mismatch at the air-human body interface and only a small fraction of the incoming energy is 
transferred into the body interior. As the result, a very high accuracy in evaluating pressure 
and velocity fields at the air-human body interface is required in order to reliably compute the 
corresponding fields inside the human interior and, in particular, in the inner ear. 

In our previous analysis [1, 2] we developed a novel, two-stage approach based on a volu¬ 
metric integral equation formulation for objects of inhomogeneous material properties. That 
enhanced approach required introducing surface distributions of the pressure and velocity fields 
associated with the external high-contrast interfaces. 

However, in the course of subsequent work, cases of slow convergence which lead to a 
reduced accuracy (e.g., for geometries containing internal high-contrast interfaces) were en¬ 
countered / identified. 

Therefore a part of our current effort was devoted to the development of a full surface 
integral equation formulation for a problem involving multiple piecewise-homogeneous material 
regions characterized by high-contrast interfaces. The constructed solver incorporates non- 
lossy, error controlled, FFT-based matrix compression which, due to its NlogN (where N is the 
number of unknowns) memory and complexity performance, allows handling large anatomically 
realistic models without compromising the accuracy. Using the developed solver we achieved 
numerically reliable results when describing wave propagation in multi-region topologies with 
individual piecewise homogeneous sub-regions. 

The developed approach/solver constitutes a significant addition with respect to both our 
previous formulation, as well as to other competitive methods. It also establishes the basis 
towards a combined surface-volume formulation which, 

- would integrate the developed volumetric and surface integral equation solvers, and, 

- would be particularly applicable to topologies consisting of adjacent, high complexity, 
strongly inhomogeneous, but relatively low contrast regions (e.g., bone, brain, ligament 
tissues) embedded in regions approximately homogeneous but characterized by high- 
contrast (e.g. skin tissue exposed to air). 

Development of such a formulation we have already initiated and to pursue its continuation 
would be of significant interest and value. 

Another ingredient critical in numerical simulations is a proper surface representation of 
the geometry of interest. Therefore significant fraction of the effort was devoted to the con¬ 
struction of a geometry model which would be anatomically faithful, would contain all the 
relevant components controlling the energy flow through the human head and, at the same 
time, would meet the stringent requirements of numerical simulations. Principal components 
of the developed human head model include: 

- the outer surface of the skin surrounding the skull and containing 

- the outer ear represented by its exterior surface, 
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- the surface of the auditory canal, 

- the tympanic membrane modeled as a finite-thickness surface, 

- the skull including the temporal bone and boundaries of the inner ear cavity: 

- the cochlea with oval and round windows, 

- the vestibule, 

- the semi-circular canals. 

- the middle ear cavity, consisting of the system of ossicles and supporting structures, 


The third main part of our effort comprised extensive numerical analysis employing the 
developed formulation and the geometry. The main results can be summarized as follows: 

- Distributions of pressure and velocity fields depend in a very nontrivial way on the 
object geometry, hence the behavior of the resulting energy flux into the inner ear may 
be difficult to predict without realistic simulations. 

- Resonances associated with the outer ear canal emanate energy to the surrounding tissues 
and provide an important contribution to the energy flux reaching the inner ear by means 
of non-airborne sound transmission. 

Comparison of the results for various head geometries shows that at relatively high 
frequencies (above about 2 kHz, corresponding the wavelengths < 15 cm) a large part 
of the non-air-borne energy transfer to the inner ear is due to the resonant behavior of 
the wave in the outer ear canals. This contribution to the energy flux can be significantly 
reduced by means of ear-plugs. However, as the frequency decreases, the contribution of 
the energy penetrating through the surfaces of the head and the skull becomes dominant; 
and this energy-transfer mechanism can only be suppressed by protecting the entire 
surface of the head. 

- The skull affects the local energy flux distribution to the inner ear, possibly through 
reflections from the interfaces between the bone and the soft tissue. 

The continuation/completion of the above described approach may constitute a solid basis 
for investigations of the effectiveness of air- and bone-conduction noise-protection devices, as 
well as the effectiveness of various hearing-enhancement aids. 

4 Publications 

• (1) ’’Formulation and applications of an integral-equation approach for solving scatter¬ 
ing problems involving an object consisting of a set of piecewise homogeneous material 
regions” 

J. Acoust. Soc. Am. Volume 130, Issue 4, pp. 2435-2435 (2011); Elizabeth Bleszynski, 
Marek Bleszynski, and Thomas Jaroszewicz 
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• (2) ’’Numerical simulation of acoustic wave interaction with inhomogeneous elastic bodies 
containing homogeneous inclusions” 

J. Acoust. Soc. Am. Volume 129, Issue 4, pp. 2382-2382 (2011); Elizabeth Bleszynski, 
Marek Bleszynski, and Thomas Jaroszewicz 

• (3) ’’Formulation and selected applications of a regularized elastodynamics integral equa¬ 
tion approach to large scale scattering problems involving inhomogeneous objects” J. 
Acoust. Soc. Am. Volume 131, Issue 4, pp. 3510-3510 (2012); Elizabeth Bleszynski, 
Marek Bleszynski, and Thomas Jaroszewicz 

• (4) ’’Reduction of Volume Integrals to Nonsingular Surface Integrals for Matrix Elements 
of Tensor and Vector Green Functions of Maxwell Equations” 

IEEE Transactions on Antennas and Propagation, vol. 61, No. 7, July 2013; Elizabeth 
Bleszynski, Marek Bleszynski, and Thomas Jaroszewicz 

• (5) “Integral-equation solver for investigation of acoustic energy flow in the human head” 

(to be submitted to JASA); Elizabeth Bleszynski, Marek Bleszynski, and Thomas Jaroszewicz 

• (6) ’’Formulation and Applications of the First and the Second Kind Elastodynamics 
Integral Equations ” 

(to be submitted to JASA); Elizabeth Bleszynski, Marek Bleszynski, and Thomas Jaroszewicz 


5 Summary of the results 

We briefly summarize here the areas of our work and our main results: 


5.1 Development of integral-equation formulation and solution algorithm 

■ We assessed merits of various forms of integral equations employing different volumetric 
and surface representations in their ability to model problems of our interest: propagation 
of an elasto-acoustic wave through the human head. 

■ We selected and constructed formulations based on first and second kind surface integral 
equations for geometries composed of multiple, homogeneous regions: 


- in acoustics, the problem is formulated for two unknown scalar fields, 


p{r ) and q(r ) 


1 dp{r ) 
p{r) dn{r ) 


where p is the pressure on an interface, p the medium density, and n the normal to 
the interface, 

- in elastodynamics, for two unknown vector fields, 


u{r ) and £(r) , 


i.e., for the displacement and traction on the interfaces. 
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■ An important element of the formulation enabling its efficient implementation was the 
reduction, by means of suitable integration by parts, of singular terms appearing in 
matrix elements arising in surface integral equations to weakly singular integrals. The 
paper describing the approach has been submitted and accepted for publication (see 

Ref. 0 ). 

■ We developed a parallel solver based on the selected formulation. The solver uses an 
FFT-based error-controlled matrix compression which 

- is suitable to problems involving highly sub-wavelength discretization (such as the 
middle and inner ear cavity problem), 

- is characterized by NlogN memory and complexity performance (where N is the 
number of unknowns), and hence allows handling large anatomically realistic models 
without compromising the accuracy. 

■ We initiated development of a combined surface-volume formulation suitable for geome¬ 
tries consisting of adjacent, high complexity, strongly inhomogeneous, but relatively low 
contrast regions (e.g., bone, brain, ligament tissues) embedded in or adjacent to regions 
approximately homogeneous but characterized by high-contrast (e.g. skin tissue exposed 
to air, lining of nasal cavity, etc.). As the first step we developed discretization methods 
for equations coupling surface and volume variables on interfaces and in inhomogeneous 
domains. 

Completion of this approach may lead to the state of the art numerical simulation tool capable 
of accurate determination of acoustic energy transfer and its deposition in geometrically very 
intricate and most vulnerable regions of a human head. 

5.2 Geometry construction 

■ Availability of a high quality geometry model of a human head of is critical importance 
in obtaining reliable results. By a high quality model we understand: 

- a model anatomically faithful, with the proper choice of relevant organs/tissues and, 
at the same time, 

- a model with geometry discretization complying with the strict requirements of the 
numerical tool. 

Therefore, geometry construction constituted a very intense part of the project. Our 
geometry model has been built by using several publicly available sources and required 
considerable work in matching the components. In its present form it contains essential 
features needed to study energy flow/transfer to the human head and to the cochlea 
region. 

■ The principal components of our human head model are: 

- the outer surface of the skin surrounding the skull and containing 
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- the outer ear represented by its exterior surface, 

- the surface of the auditory canal, and 

- the tympanic membrane modeled as a finite-thickness surface, 

- the skull, described by surfaces of the bones and containing 

- the required boundaries of the middle ear cavity, 

- the boundaries of the inner ear cavity (the cochlea, the vestibule, and the semi¬ 
circular canals), 

- the middle ear cavity, consisting of the system of ossicles and supporting structures, 

■ We stress that there are several imperative requirement in constructing a geometry repre¬ 
sentation suitable for numerical simulations which make it a tedious and time consuming 
process. They include: 

- nearly uniform discretization with local variations reproducing small geometry de¬ 
tails, 

- good quality (aspect ratio) triangular mesh, 

- closed surfaces with no intersecting or overlapping patches and, 

- grid connectivity on joints between patches. 

Below we show some representative examples of the head geometry used in our simula¬ 
tions. 

Please zoom the pdf file when you are viewing the figures below to see better the geometry 
details and the mesh quality. 


6 



Figure 1: The discretized skin geometry. 



Figure 2: A detail of the discretized skin geometry viewed from the inside - the outer meatus. 
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Figure 3: A view of the skull geometry modeled as a closed surface with the ear canal and the 
ossicles in the middle ear cavity (the eardrum was removed for a better view). 



Figure 4: A close view at the outer ear canal, the middle ear cavity and the ossicles (the 
eardrum was removed for a better view). 
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Figure 5: A view of the skull geometry modeled as a closed surface from the inside - note the 
temporal bone and the inner meatus. 



Figure 6: A detail of the skull surface: the inner meatus, cochlea, semicircular canals, the 
middle ear cavity). The triangulation used for this particular component of the geometry 
resulted in about 280, 000 unknowns. 
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Figure 7: Another detail of the skull surface: the middle ear cavity with the attic, the semi¬ 
circular canals, and the ossicles (with the eardrum removed for a better view). 
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Figure 8: A part of the discretized head geometry model in the vicinity of the left ear, shown as 
a coronal section seen from the back. The dark surfaces visible in the cross-section are interior 
sides of the boundaries of the bone and the outer boundary of the head (the skin). The space 
between the skin and skull, and the interior of the skull are filled with soft tissue. Since the 
clipping plane intersects the outer ear canal, its interior is partly exposed. 
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■ Although the degree of detail and complexity of our model allows us to carry out rea¬ 
sonably realistic numerical simulations, it would be beneficial to further enhance the 
geometry model by including, in particular, such additional elements as 

- the nasal cavity, 

- the brain tissue, the latter characterized by nontrivial elasto-viscous properties. 

Such components should be carefully and in as much as possible complete manner taken 
into account, as they may play a significant role in controlling the distribution of the 
energy flow into the inner ear and therefore further improve fidelity of the simulations. 

An important and very valuable aspect of the created geometry is that, due to its perfect 
connectivity, the shapes and sizes of individual organs can be locally deformed (rescaled) leading 
to different head geometry realizations. Hence, sensitivity of the energy flow and energy flux 
distribution on the geometrical details of the human head can be studied. 

5.3 Numerical simulations and assessment of the relative importance of 
the bone conduction mechanism 

5.3.1 Overview 

As we mentioned in the Introduction, the acoustic and elastic wave propagation problems we 
are investigating involve penetration of waves from a low-density medium (air) to the interior 
of an object consisting of high-density tissues. 

This category of problems have been relatively little explored in the past, either theoreti¬ 
cally or in applications. Most of the theoretical work concentrated on essentially single-region 
problems involving hard- or soft-surface or impedance boundary conditions, while the transmis¬ 
sion boundary problems have been receiving much less attention. Similarly, typical applications 
were either (1) dealing with much lower contrasts than in our case (for instance, in marine 
applications, rock- or metal-to-water density ratios not exceeding about 10, compared to our 
ratios > 1000), or (2) considering wave sources embedded in the dense medium rather than in 
the surrounding air, as is the case in seismic or vibration mechanics problems. 

In view of this status of the current research, our problem required a considerable effort to 
develop an insight into the mechanisms of wave penetration through high-contrast interfaces 
and its propagation in the dense medium. In particular, our extensive numerical simulations 
allowed us to to establish a general picture of energy flow inside the human head model and 
identify the main mechanisms controlling its distribution (including reflection from and trans¬ 
mission through high contrast interfaces, as well as resonances in air cavities); we were also 
able to obtain quantitative estimates of the amount of energy (relative to the incident energy 
flux density) reaching the inner ear by means of bone- and soft-tissue conduction. We stress 
here that, although the overall level of the energy flux entering the cochlear cavity is rela¬ 
tively insensitive to the geometry and material parameters variations, its detailed distribution 
and even orientation do depend on the geometry in a nontrivial way and would have been 
impossible to predict without the actual computations. 
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The main results of the simulations are presented below. Additional results are contained 
in the attached draft of a paper (Appendix B). 

5.3.2 Energy flux as a measure of the tissue conduction effects 

As a quantitative measure of the amount of bone-conducted acoustic energy reaching the inner 
ear we compute the integral of the absolute value of the average energy flux density F over the 
surface S of the cochlear cavity, 




(5.1) 


where \S\ is the area of the surface. 

5.3.3 Bone-conducted energy flow into inner ear in the presence and absence of 
the external ear canal 

We present here results for the energy flow distribution in two models, both including the full 
skull and middle/inner ear cavity, but differing in the presence or absence of the outer ear 
canal. The latter model can be considered a realization of a perfect ear-plug , with which the 
outer auditory canal is effectively closed. 

In both models we assumes typical approximate values of the relative densities and the 
refractive indices, p/p 0 = 1000 and n = 0.2 for the soft tissues, and p/p 0 = 2000 and n = 0.4 
for the bone. The model is being subject to a pressure plane wave of unit amplitude, incident 
from the left side of the head. In most of the examples we carried out computations in a wide 
range of wavelengths, from A = 3cm (about 10kHz) to A = 80cm (about 400Hz); A = 27i/k 0 
is always the wavelengths in the air. 

We note that in these simulations we excluded the impedance-matching mechanism of 
the middle ear; hence, the entire energy arriving at the inner ear has propagated exclusively 
through the bone and other tissues, and not through the air-conduction pathway. 

Distribution of the energy flux in the temporal bone area for A = 5 cm is visualized in 
Fig. 9(a). The flux density has large values on the walls of the outer ear canal and the 
energy evidently emanates from there to the surrounding tissues. In fact, by comparing this 
distribution with that for the alternative model we find that the presence of the outer ear canal 
increases the energy flux through the inner ear by about a factor of 10. 

An interesting feature is that in different parts of the canal the energy flows from the air 
into the tissue or in the opposite direction. Moreover, we found that the energy flow changes 
direction several times with the changing wavelength in the considered range. Such a behavior 
is characteristic of a standing acoustic wave being formed in the ear canal [3]. 

A more systematic and quantitative comparison of the models is provided in Fig. 9(b) 
which shows the average relative energy flux through the cochlea (Eq. (5.1)) as a function of 
the wavelength of the incident wave. At larger wavelengths the presence of the outer auditory 
canal increases the energy flux by about a factor of two. At smaller wavelengths the effect is 
larger and exhibits a resonant behavior, manifesting itself by the two peaks at A ^ 5 cm and 
A ~ 15 cm. 
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(b) 


Figure 9: (a) Distribution of the energy flux density on the material interfaces in the vicinity of 
the left temporal bone, for A = 5 cm. Positive (red) and negative (blue) flux indicates energy 
flowing into and out of the bone, (b) The average energy flux density through the surface of 
the cochlear cavity (relative to the incident wave flux density), as a function of the wavelength, 
for the model with the skull and soft tissues. 


5.3.4 The role of acoustic bone conduction vs. soft tissue conduction 

In order to assess the significance of bone conduction vs. soft-tissue conduction, we also com¬ 
puted the amount of energy flowing through the cochlea in a model of the human head filled 
entirely with a homogeneous material, without the skull structure. In this comparison we con¬ 
sidered, as before, two models of the outer head surface: without outer ear canals (i.e., with 
“perfect ear-plugs”) and with ear canals. 

We first consider the model without the outer ear canals, in which the entire head is filled 
with the soft tissue (p/p 0 = 1000, n = 0.2). We solve the transmission problems with the 
outer head surface (the “skin”) and a part of the middle- and inner-ear ear cavity, including 
the cochlear cavity. The presence of the latter surface does not affect the solution itself, but 
allows us to compute the average energy flux 

* 8 :=±J#r\F(r)\, := Kp(r) Im h( r ) g nP*( r )} > ( 5 - 2 ) 

through the same part S of the cochlea cavity surface as used before. (Fig. 10). 
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(a) 

Figure 10: (a) A part of the discretized head geometry model in the vicinity of the left ear, 
shown as a coronal section seen from the back. The dark surfaces visible in the cross-section 
are interior sides of the boundaries of the bone and the outer boundary of the head (the skin). 
The space between the skin and skull, and the interior of the skull are filled with soft tissue. 
Since the clipping plane intersects the outer ear canal, its interior is partly exposed, (b) Part of 
the discretized geometry containing the left middle- and inner-ear cavity, viewed from below. 
The surface of the cochlear cavity is marked with S. 

The results are shown in Fig. 11. It can be seen that the reduction of the flux due to the 
presence of the skull is not large, of the order of 40%, and can be plausibly attributed to an 
additional reflection of the sound wave from the interfaces between the skull and soft tissues. 

Results for the analogous model with the outer ear canals are shown in Fig. 11(11). Again, 
the flux densities in the the three models are not strongly affected by the presence of the skull. 
In particular, the persistent resonant structure suggests that it is due rather to the shape of 
the head itself than to the presence of the skull. 
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Figure 11: The average relative energy flux density through the the cochlear cavity as a function 
of the wavelength, for the model without the skull (b), compared to the result for the reference 
model (a) with the skull. Plots (I) and (II) correspond to head surface models without and 
with the outer ear canals. 


The above comparisons of the absolute energy flux densities might seem to suggest the skull 
does not play an important role in sound conduction. However, a more careful examination of 
the spatial distribution and orientation of the energy flow shows this is not the case: 

As an example of such an analysis, we plot in Figs. 12 flux distributions on the middle- and 
inner-ear cavity at the wavelength 12cm (in the resonance region), for the model with (a) and 
without (b) the skull, but in the presence of the outer ear canals. We find that the direction of 
the energy flow in these models is almost exactly reversed . This phenomenon (which is, in fact, 
observed for most of the considered wavelengths) has to be attributed to interactions (perhaps 
reflections) involving the skull bones. 
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Figure 12: Energy flux density distributions on the left middle- and inner-ear cavity for wave¬ 
length 12 cm computed for the full model (a) and for the model without the skull (b). The 
cavity surface is viewed from the top; the canal directed to the left and up is the inner meatus. 
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5.3.5 Summary of simulation results 

In summary: 

- We were able to establish a general picture of energy flow inside the head model and 
identify the main non-air conduction mechanisms (in particular, transmission through 
and reflection from high-contrast interfaces, as well as resonances in air cavities). 

- We found that the overall magnitude of the energy flux reaching the inner ear is relatively 
insensitive (within a factor x to the details of the geometry and even the presence or 
absence of the skull (hence, sound conduction through soft tissues is also significant). 

- However, the more detailed energy flux distribution and even the flux orientation proved 
to be highly sensitive to the shape of the interfaces and to the difference of material 
properties between the bone and the soft tissues. For instance, the presence of the skull 
bones may reverse the direction of energy flux near the middle ear. Such effects may 
be of importance in a more detailed analysis of the physiological effects of the energy 
delivered to the inner ear. 

- We established the frequency dependence of the energy reaching the inner ear through 
the tissues. We found that at higher frequencies (few centimeter wavelengths in air) a 
large part of that energy emanates from the outer ear canal, which supports a resonant 
standing wave. However, at lower frequencies, the energy arriving at the inner ear is 
dominated by transmission through a large area of the head surface. 

- The last result suggests that the higher frequency noise components may by significantly 
reduced by blocking the ear canals, but lower frequencies can only be suppressed by 
protecting the entire surface of the head. 


6 Open problems requiring possible future research 

Although we believe we achieved a substantial degree of understanding of the sound propaga¬ 
tion mechanisms in the human head, we also encountered challenges requiring further investi¬ 
gation and offering, potentially, interesting and useful results. We list some of them: 

- The role of large cavities in the head (in particular, the nasal cavity) should be inves¬ 
tigated more thoroughly. Simulations should include possible shape variations (hence 
possibly different resonant behavior) and should establish the significance of the cavities’ 
communication with the outer atmosphere. 

- The skull bones are known to contain a multitude of small air cavities effectively forming 
a porous bone structure; their presence may significantly affect wave scattering. Because 
of their complexity and the lack of a reliable geometry model, the small cavities were not 
included in our simulations. It might be possible to take them into account by means of 
evaluating effective medium parameters, which would then probably exhibit substantial 
dispersion. 
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The brain tissue is known to exhibit significant viscosity. Its effects in energy dissipation 
should be investigated. 

We found that although the overall magnitude of the tissue-conducted energy reaching 
the inner ear is relatively stable, the details of its spatial distribution strongly depend on 
the head geometry. It would be of interest to investigate more thoroughly these effects 
and their possible physiological role. 

In our simulations we found that high-contrast material layers (such as skull bones in 
air) reflect most of the acoustic energy and transmit very little, even if they are thin 
compared to the wavelength. This behavior suggests possible efficient noise protection 
devices based on high-contrast layered structures. A detailed investigation conducted 
with the developed simulation tools may proof to be of significant interest and usefulness. 


17 


Appendices 

A Details of the formulation 

We provide here the main elements of the theory and some details of implementation of our 
solver for elastodynamics problems. More information on the formulation of the acoustic 
integral-equations and their implementation in our current solver are given in the draft of a 
paper contained in Appendix B. 


A.l Formulation of elastodynamics integral equations and algorithm devel¬ 
opment 

In this Section we provide details on the elasto-dynamics integral equations for multi-domain 
problems we constructed and partially implemented. We start with the Lame equation for 
the displacement in elastic medium subject to an external plane wave excitation F(r) The 
displacement u{r) satisfies the equations 

{(A + /i) V r ®V r + /iVj + pu 2 }u(r) = F(r), (A.l) 

where A and p are the (position dependent) Lame parameters, p is the medium density and uo 
is the incident wave frequency. 

The homogeneous medium Green function of the Lame equation is the second rank sym¬ 
metric tensor 


G{R) = - g s {R) I + -L ® Vr [g s (R) - gc(R )] , 




gk 2 s 


(A.2) 


where 


gc{R) = 
gs(R ) = 


Q i kcR 

4tt R ’ 

a i ksR 


4t tR ’ 

R = x — y . 


(A.3) 


The two wave-numbers, 


k c — 


UJ 


k s — 


uj 


(A.4) 


are related to the longitudinal (compressional) and transverse (shear) wave speeds, 


c C ~ 


' A -j- 2p 


P 


(A.5) 


We consider a single material region Q bounded by dQ. By assuming that, inside 17, the 
displacement u(x) satisfies the Lame equation (C.14), we may construct the following integral 
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representation for the displacement field , which allows us to find its value at any point 
x in Q in terms of an integral of two fields, u(x) and i(x), on the region boundary: 


u(x) = d 2 y[u(y) ■ T(x,y)+t(y) ■ G{x - y)] . (A.6) 

Jan 

The symbols appearing in the integrand of the above integral representation are as follows: 

• u(y) is the displacement vector field, 

• t(y) is the traction vector field related to the stress tensor r(y) as follows: 

t(y) = n(y) ■ f(y) , 

f(y) = A 7 V y • u(y) + y [V y <g> u(y) + u(y) <g> V y ] , 

t(y) = A n(y) [V y • u(y)] + y{[n(y) ■ V„] u(y) + V y [n(y) ■ u(y)}} , (A.7) 

Tij{y) — Cijkm d k u m (y) , 

C’ijkl = ^ A // (8j k Sji + 8j[ 8j k j , 

• G(x — y) is a second rank symmetric tensor, the Green function of the Lame equation 
given by (C.2), 

• r(x,y) is a second rank non-symmetric tensor, related to the third rank stress tensor 
Green function U(x) (symmetric in its first two indices) 

r{x,y) = n(y) • U(R) , 
d 

rij(x,y) — Gip(x y ) CjkpqTik^y) { 71 ( 3 /) • S{x y)}ij ? 

^Uq 

Rjk{x y) = ni(y ) Cij m i d m Gi k (R) , 

= A <% c> m (R) + y [di G jk (R) + dj G ik (.R)] = c> m G /fc (iZ) , ( ' 8) 
y) = My) Eijk(R) 1 
@jk(x,y) = E ijk (R)rii(x) , 

E(R) = - E(-R) . 

The tensor Green functions G(R) and E(R) satisfy the following equations: 

[(A + y) S7 R ®V R + yV 2 R + P uj 2 }G{R) = - S(R) I , 

V fi • E(R) + pu 2 G(R) = - 8(R) I , 

E(R) • V fi + pa? G(R) - A (V fl • -V R ®)V R • G(R) = - 8(R) / , 

di [C ijkl d k Gi m (R)\ + pa? Gj m (R) = - 8 jm S(R) , ( A - 9 ) 

di [Eij m (R)\ + pa?Gj m (R ) = — 8j m 8(R) , 

• G(R) = \ 2 V r9c (R) , 

paj z 
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while t(x) satisfies the equation 


V r • r(r) + poj 2 u(r ) = 0 , 
i[V r <g> u(r) + u(r ) (g> V r ] - D ijld m{r) = 0 , 


(A.10) 


with 



(A.ll) 


is the fourth-rank elastic stiffness tensor and 


Dijkl (^ik 3ji T fiil fijk) 2 ^ ^ij^kl 


(A.12) 


An additional useful relation is 


Cijki &k di Gj m (yR) Cijhi d{ Uj(x ) dfo G/ m (-R) — 0 . 


(A.13) 


A.2 Construction of coupled surface integral equations for displacement and 

surface-traction fields in piecewise homogeneous media 

Surface integral equations or boundary integral equations (BIEs) are applicable to piecewise 
homogeneous materials, and provide solutions for the displacement and traction fields defined 
on interfaces 5 mri separating different material regions (Fig. 13). One of these regions, 
is the unbounded background medium (air). Fields in the individual regions are described in 
terms of the appropriate Green functions for elastic materials. The displacement and traction 
fields are assumed to be continuous across the interfaces (i.e., to satisfy transmission boundary 
conditions). 


,/ 


r 



Figure 13: A schematic representation of regions Q and interfaces S appearing in surface 
integral equations . 

As an example, we give below explicit expressions for two alternative systems of surface 
integral equations describing for the above transmission problem. Such systems can be then 
obtained in the following three steps: 
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(a) using a conventional representation theorem (found, e.g., in Ref. [4]) for each material 
region i? for the displacement field in this region written in terms of the surface integrals 
over dfl with displacement and traction fields, u(r) and t(r) = n(r) • f (r), playing the 
role of surface sources, 

u(x) = [ d 2 y[r^(x,y)-u(y)+ G^(x,y)-t(y)] ; (A.14) 

JdQ 

(b) imposing continuity boundary conditions on the u(r) and t(r) values on each interface 
(oriented surface) S m n separating the regions Q m (on the negative side of the interface) 
and fl n (on its positive side), 


u m {r) = u n (r) , (A.15a) 

tm(r) = t n (r) ; (A.15b) 


and 

(c) writing two suitable equations following from the boundary conditions for each interface 

Smn- 

The choice of the integral equations is not unique. It is possible to form several different 
sets of integral equations by taking different linear combinations of equations representing the 
boundary conditions (A. 15). While all such integral equations are theoretically equivalent, they 
tend to differ in terms of accuracy, computational resources needed, and solution convergence. 
We considered two representative choices marked below as (i) and (ii) and, after examining 
them, we decided to choose the second, non-conventional, integral equation set. 

The set (i) of integral equations are obtained by imposing, on each interface S mn , the 
transmission conditions on the displacements only (Eq. (A. 15b)), but expressing them by means 
of the two representation formulae, Eq. (A.14) and its counterpart for the other region, 
adjacent to the interface. The resulting system of six equations for six unknown components 
of u(x) and t(x) is 

\ «(*) + [ d 2 y [r£(x, y) ■ u(y) + G^(x, y) ■ t(y)] 

J 8 mn 

- E / d2y [ r m( x i y ) • u (y ) + G m ( x > y) • %)] = <L 0 u m (x) for * e s mn (i-a), 

S. £df 2 Sim 

im m 

(A.16a) 

5«(*)- [ d 2 y[r^(x,y)-u(y) + Gl(x,y)-t(y)] 

Js 

mn 

+ E / d2y i r n( X ’ y) ■ u (y ) + G n (*, y) ■ t(y)\ = Ko uin ( X ) for * e S mn ( i_b ) • 

S .edf 2 S n j 
nj n 
j^m 

(A. 16b) 
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With reference to Fig. 13, the first of the above integral equations represents contributions to 
the displacement field u on the interface S mn due to the displacement and traction fields u and 
t on the same interface (the first integral) and on remaining interfaces, S^ m , forming boundaries 
of the region f2 m with other regions i ^ n; and similarly for the second equation, which 
represents contributions of the boundaries of the region fi n . The r.h.s.s of the above equations 
are the incident fields due to distant sources in the region hence the delta-functions 5 m o 
and 5 n0 . 

We obtain the second, alternative, set (ii) of the integral equations by imposing the bound¬ 
ary conditions (A. 15) on both the displacement and traction fields. We again express the 
displacements u by means of the representation formula (A. 14) and its Q n counterpart. For 
the traction fields t we use the same representation formulae to which we apply the trac¬ 
tion field operator n(x) • r[u(x)\. The resulting system (ii) of six equations for six unknown 
components of u(x) and t(x) is now 

[ d2 y {l r m(x,y) + (x,y)\ ■ u(y) + [G^(x,y) + G^(x,y)] •%)} 

J s mn 

~ J2 f d 2 y{r£(x,y)-u(y) + Gl l (x,y)-t(y)} 

S. edn S im 

im m 

iy^n 

~ J d 2 y{ r „(x,y)}-u(y) + Gl(x,y)-t(y)} 

S .ed£2 ^nj 
nj n 

j^m 

= l 5 mO ~ 5 no] «“(*) fOT * 6 S mn (M > 

(A.17a) 

/ d 2 y {[W^(X; y ) + Wj(x, y)] ■ u(y) + [$ T (x, y) + ( Pj L (. x , y)} ■ t(y )} 

JS 

mn 

~ J2 f d 2 y{w£(x,y)-u(y) + $l l (x,y)-t(y)} 

S. EdO S im 

im m 

iy^n 

~ j d 2 y{W^(x,y)]-u(y)+$l(x,y)-t(y)} 

- ^nj 

= i 5 mO ~ 5 no] *“(*) fOT * 6 S mn (^0 > 


S .edf2 nj 
nj n 

j^m 


where, in the second equation, the boundary-operator kernels are defined by 


(A. 17b) 


W (*, y) ■ u{y) = n(x) ■ t[ r(x, y) ■ u(y)} , 

®(x, y) • t(y) = n(x) ■ t[G(x , y) ■ t(y)\ . 

As before, the operators with indices m and n describe propagation of the displacement and 
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stress fields in the regions fi m and fi n . The details of the construction and regularization of the 
kernels G^, G, Wi, and ^ are described in a paper we are preparing for publication in JASA. 

The main advantage of the set (ii) of the integral equations is that it results in higher 
solution accuracy at material discontinuities; it is achieved through explicit imposition of 
boundary conditions on the displacement and the traction fields, rather than by means of the 
weaker integral relations used in the formulation (i). We also find that the second choice, 
(ii), yields a better solution convergence in high contrast problems. Theoretical aspects of the 
differences between the two above-mentioned choices, (i) and (ii), will be more fully investigated 
in the near future. 

A.3 Matrix compression 

The underlying element of our approach is the Fast Fourier Transform (FFT)-based AIM 
matrix compression method, initially developed in the context of electromagnetics for solving 
large-scale problems, and described in detail in Ref. [5]. Adaptation of this formulation to 
large-scale acoustic problems [1, 2] was the initial step of our effort. 

The main reason for choosing the FFT-based compression method, rather than other com¬ 
pression techniques, is that it provides superior efficiency in the treatment of both volumetric 
and surface integral equations, particularly for sub-wavelength problems (geometrical elements 
used in objects’ discretizations - triangular facets and tetrahedrons - are much smaller than 
the wavelength). We note that sub-wavelength geometry regions constitute dominant portion 
of anatomically realistic head geometry models. 
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B Draft of the paper 

“Integral-equation solver for investigation of acoustic energy 
flow in the human head” 

E. Bleszynski, M. Bleszynski, and T. Jaroszewicz 
Monopole Research, Thousand Oaks, CA 91360 
Abstract 

We describe selected aspects of the development and applications of integral equation 
based solver designed for large scale numerical simulation of interaction of sound waves with 
complex heterogeneous media. A particular application we are concerned with is modeling 
of the propagation of the waves inside the human head subject to a pressure wave in the 
surrounding air; the goal of this investigation is to estimate the amount of energy trans¬ 
ferred to the inner ear by air- and bone-conduction mechanisms. The considered problem 
is characterized by a high density contrast (biological tissues vs. air), which causes the 
dominant part of the incident energy to be reflected, and makes the computation of the 
penetrating energy flux rather nontrivial. In our previous work we developed and employed 
a computational model based on volumetric integral-equations of acoustic/elastodynamics, 
with the object modeled by means of a tetrahedral mesh. Since in this formulation the 
equations were only providing the pressure (or displacement) distribution, computation of 
the energy flux required numerical (finite-difference) evaluation of the normal derivative of 
the pressure (or the traction), which we did not find sufficiently reliable, especially in the 
geometrically complex region of the model near the inner ear. The present formulation is 
based on surface integral equations, which provide simultaneously solutions for the pressure 
and its normal derivative (or for the displacement and the traction). The equations are 
discretized in terms of piecewise-constant and piecewise-linear basis functions supported 
on facets of a triangulated surface model representing the human head including the skull, 
the soft tissues, and the geometrical details of the middle and inner ear. The resulting stiff¬ 
ness matrix is compressed by means of the Fast Fourier Transform (FFT)-based algorithm, 
implemented in distributed-memory parallel processing environment, which allows han¬ 
dling of problems involving several million unknowns. Particular attention was paid to the 
appropriate choice of the compression parameters in the practically available anatomical 
models characterized by large differences of the mesh densities. We also present results em¬ 
ploying several alternative integral-equation formulations of the multi-region transmission 
problem, giving rise to either first- or second-kind equations, with different conditioning 
properties. In this context we analyzed some pertinent preconditioning techniques, includ¬ 
ing a Calderon-type preconditioner for the first-kind and an additive-type preconditioner 
for the second-kind equations in the high-contrast problems. 

B.l Introduction 

The purpose of this work is to analyze some aspects of wave propagation and the resulting 
energy flow in a human head subject to an incident acoustic wave propagating in the surround¬ 
ing air. This analysis is carried with the goal of estimating the amounts of energy reaching 
the inner ear, due to various sound wave conduction mechanisms. In the following we briefly 
discuss difficulties arising in this investigation and the solutions implemented in our work. 
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Mechanisms of energy transport in the human head. It is known that the normal 
airborne mechanism of acoustic energy transfer to the cochlea - through the outer auditory 
canal, the tympanic membrane, and the middle-ear system of ossicles - is not the only one 
in operation. Alternative energy flow paths, in particular due to bone conduction, are being 
utilized in, e.g., in hearing-enhancement devices. The same mechanisms, however, can also 
transfer to the inner ear the energy due to an intense noise; since those pathways circumvent 
the usual airborne sound transmission channel, they may potentially cause inner ear damage 
in spite of application of conventional protective devices, such as ear plugs. 

Significance of high medium contrasts in evaluation of energy transfer. A charac¬ 
teristic feature of the considered problem of acoustic wave propagation is a very high density 
contrast (density ratio of the order of 1000) between the biological tissues and the surrounding 
air; at the same time, the refractive indices of those media differ by only a modest factor 
of at most 5. The high density contrast causes a large impedance mismatch and greatly re¬ 
duces the energy transfer through the air-tissue interface. This circumstance is the reason why 
the ear structure of animals living in the atmosphere evolved to form an exceedingly sensi¬ 
tive impedance-matching mechanism consisting of the tympanic membrane and the middle-ear 
system. 

In this context, we note that the area of acoustics and elasticity we are concerned with - 
high-density objects subject to waves incident from the surrounding low-density medium - has 
been relatively little explored. For instance, 

- Marine applications involve interaction of waves in water with denser materials (metallic 
objects, rock), but the density contrast is much lower than in our case; e.g., steel-to-water 
density ratio is about 8. 

- In geophysical (seismic) applications the object of interest (the Earth) is surrounded by 
air, but wave propagation through the atmosphere is not being considered at all. 

- Mechanical engineering problems often involve interaction of vibrating structures with 
air; however, it is the vibrating structure (e.g., a body of a vehicle) which is the source 
of sound waves, and the effect of the acoustic waves on that structure is negligible. 

- Finally, relatively little theoretical work has been done on transmission problems involv¬ 
ing high-contrast media. 

Computational aspects. In view of the high density contrasts and the the resulting impedance 
mismatch, reliable computation of the amount of energy reaching the inner ear becomes a 
highly demanding task. This energy is relatively small and its value depends on the energy 
flux though the high-contrast air-tissue interface (including effects of many air cavities present 
in the head). The energy flux, in turn, is proportional to the product of two physical quan¬ 
tities on the outer side of the interface: the pressure p and its normal derivative dp/dn (the 
latter quantity is, essentially, the normal component of the velocity of the medium). Because 
of the large density contrast, the problem of the tissue-air interface is nearly that of the hard 
surface: the value of p on the surface is roughly of the same magnitude as the pressure in 
the surrounding air, but the normal derivative dp/dn is much smaller than in the surrounding 
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space, and its evaluation is not entirely straightforward. In addition, dp/dn has to be known 
to a sufficient accuracy, since the relative energy flux value is linear in that quantity and the 
same is true for the error in the energy flux. 

Past and present approaches. In our previous analysis [1, 2] we developed and imple¬ 
mented an approach based on a volumetric integral equation for an inhomogeneous material. 
Because of the high-contrast nature of the considered problems, that approach required, ef¬ 
fectively, introducing surface distributions of the pressure and velocity fields associated with 
high-contrast interfaces. 

In this work we extended the surface-related elements of our previous approach to a full 
surface-equation formulation for a problem involving multiple piecewise-homogeneous mate¬ 
rial regions and high-contrast interfaces. An additional motivation for this development was 
availability of detailed surface geometry data for the middle and inner ear structures: While 
moderate-resolution voxel-type data are readily accessible for entire human bodies (e.g., [6]), 
this is not the case for the more detailed skull and inner ear regions. Furthermore, a practical 
obstacle in utilizing voxel data is the necessity of “segmentation”, i.e., identification of tissue 
types - a difficult and time-consuming process that can be only partly automated. At the same 
time, more precise middle- and inner-ear data based on micro-magnetic-resonance imaging are 
typically available in the form of reconstructed interfaces between tissues. 

Further, having implemented both volume and surface solution methods, we now envisage 
a comprehensive approach combining them into a full solvers based on a system of mutually 
coupled equations. It could be, in particular, applied to models consisting of both volumetric 
regions (with a relatively “coarse” discretization of one or few millimeters) and of detailed 
interfaces discretized with sub-millimeter resolution. In fact, in the surface geometries we have 
been using in the computations reported here, triangle sizes span two order of magnitude, from 
about 5 mm to 0.05 mm = 50 /mi. 

Contents of the paper. The paper is organized as follows: In Section B.2 we present an 
integral-equation formulation applicable to two types of multi-domain problems which involve 
either a collection of piecewise homogeneous material domains or a collection of both piece- 
wise homogeneous and inhomogeneous domains. In Section B.3 we discuss some aspects of 
discretization of the equations, in particular in the presence of junctions of several interfaces. 
We also describe some basic elements of the FFT-based matrix compression and matrix-vector 
product acceleration, which allow us to solve large scale problems. In Section B.5 we present 
results of representative numerical simulations for a geometry consisting of a human head con¬ 
taining a detailed geometry representation of the outer, middle, and inner ear. Results of such 
simulations may provide an insight into physics of sound wave propagation and energy flow 
to cochlea cavity obtained with nontrivial geometries of a realistic human head model and 
demonstrate the present capabilities of our solver. 

B.2 Integral equations for the multi-region transmission problem in acous¬ 
tics 

Simulating sound propagation in an anatomically realistic model of the human head requires 
modeling of topologically complex structure of material regions and interfaces. The driving 
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physical mechanism controlling the energy flow through the model is the variation of densi¬ 
ties between the material regions (modeled as homogeneous domains). Therefore, in order 
to obtain a reliable prediction for the energy transfer across the interfaces it is necessary to 
accurately discretize the pressure and its normal derivative on complicated surface geometries. 
The surfaces may be nested inside one another, a number of surfaces may be enclosed in an¬ 
other one, and several surfaces may form a “junction”, i.e., share a common boundary line or 
a boundary point (such structures are sometimes referred to as “triple-points” or “multiple 
points” [7]). Some of that complexity of the set of surfaces may also be introduced for con¬ 
venience of geometry modeling (e.g., a complex surface may be built of a number of simpler 
patches). 

A representative geometry structure of our interest is depicted schematically in Fig. 56. It 
involves a number of soft tissues, brain, and the skull with air cavities. The interfaces between 
various materials form a number of junctions, which, in this case, are common boundaries of 
triplets of interfaces. 



junctions 

Figure 14: A simplified schematics of the topological structure of regions and interfaces in head 
model. 


Since we are primarily interested in computing the energy flow, it is convenient to use as the 
unknowns the fields appearing in the expression for the energy flux (Eq. (B.4)): the pressure 
p and its normal derivative d n p (which we call, for simplicity, the velocity field). We obtain 
the surface integral equations in the following steps: 

- We write the representation formulae for the pressure fields p on boundaries of material 
region i? m , in terms of the values of p and d n p on the same region boundaries. 

- Write the representation formulae for the velocity field d n p of the boundaries df2 m , again 
in terms of the boundary values of p and d n p. This representation formula is obtained 
by taking the normal derivative of the previous formula. 

- We impose boundary conditions (Eq. (B.l)) on the boundaries (i.e., region-region inter¬ 
faces). These conditions relate both field p and d n p being represented by the integrals 
of the representation formulae and the quantities appearing in the integrands. 


27 




- We form appropriate linear combinations of the representation formulae. Here we may 
choose to use representation formulae for both p and d n p), and obtain in this way 
first-kind (elliptic) equations; or we may utilize only the representation formulae for the 
pressure and obtain second-kind equations. In both cases the unknowns in the equations 
are the quantities appearing in the integrands of the representation formulae, i.e., the 
boundary values of p and d n p. 


B.2.1 Differential equations of acoustics and boundary conditions 

In this section we briefly summarize here the differential equations governing acoustic fields, 
which will be used to obtain the integral equations described in Section B.2. We consider 
a bounded spatial domain ft filled with a material characterized by (generally) position- 
dependent density and compressibility p(r) and n{r). The surrounding unbounded “outer 
region” i2 0 is assumed to be homogeneous and described by the “background” parameters p 0 
and ft 0 . Correspondingly, in a time-harmonic problem with frequency uo the wave number in 
the background medium is k 0 := yjp 0 uo . The Euler equation for the pressure field p(r) 
(known, in this form, as the Bergmann equation [8, 4]) is 

= (B.i) 

Po \p(r) J p 0 K 0 

Equation (B.I) for an inhomogeneous object can be solved as the volumetric Lippmann- 
Schwinger integral equation. We followed this venue in our previous work [1, 2], where we 
developed efficient numerical techniques for handling of situations relevant in the considered 
applications, i.e., discontinuous material parameters [9] and high-contrast interfaces. 

The general result following from our previous investigations was that the energy flow 
through the object was essentially smooth even in inhomogeneous regions, as long as the 
material density was varying in a range typical of biological tissues (about a factor of two); 
however, the energy flux through high-contrast interfaces was strongly dependent on the shape 
and geometrical details of the interface. 

Since the behavior of the energy flux is of main interest in the present application, we 
concentrate here on a surface form of the integral equations, assuming the object is modeled as 
a piecewise-homogeneous structure, isuch that in each homogeneous sub-region f2 m described 
by the parameters (p m , n m ) the pressure field p m satisfies the Helmholtz equation 

(V 2 + k m ) pir) = 0 for r e . (B.2) 


with the wave number k m := yjp m n rn cj . The original equation (B.I) implies then the trans¬ 
mission boundary conditions on the interfaces of regions, 


p{r) 


and 


Po dp{r) 
p(r ) dn(r ) 


~T\ n (r) ■ Vp(r) 
p(r) 


continuous across the interface ; (B.3) 


here and below n is the unit normal to the interface, and normal derivatives are denoted with 
d r ,. 
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B.2.2 Surface integral equations for p(r) and dp(r)/dn(r) for a single material 
body 

In this Section we present a brief derivation of the surface integral equations for a problem 
involving a single material region. The two unknown functions appearing in such integral 
equations are chosen to be as pressure p(r) and the normal derivative of pressure defined on 
material interfaces. 

We recount here the known basic forms of the surface integral equations for the transmission 
problem involving only a single body occupying a domain Q m (Fig. 15). 



Figure 15: 
ary T m is 
“negative” 


The simplest transmission problem with a single material region i? m . The bound- 
oriented according to the direction of the normal, and are its “positive” and 
sides. 


While there is some freedom in choosing the unknown quantities in integral-equation formu¬ 
lations, we opted for equations obtained by means of a direct method , in which the unknowns 
are the values of the physical quantities - the pressure and its normal derivative (i.e., velocity) 
on the surfaces. An advantage of such a form of integral equations is that the unknown fields 
appear directly in the expression for the density of the energy flux through a surface, 

F(r) := - Im {p*(r) d n p(r)} ; (B.4) 

this quantity is normalized to be unity for a unit-amplitude plane pressure wave incident 
normally on the surface, and, in view of the boundary conditions (B.2), is continuous across 
the interfaces. 

It is convenient to use in this formulation not directly the normal derivatives of the pressure 
field, but rather those derivatives divided by the density of the medium. More precisely, for a 
region f2 m and its boundary we to define the boundary fields 


q(r) := 


Po_ dpjr ) 
Pm dn(r) 


(B.5) 


for r approaching the boundary from inside F? m . According to the boundary conditions (B.3), 
those fields are simply continuous across the interfaces. As a result, the energy flux density 
(B.4) through any interface can be expressed in terms of fields p and q as 


F(r) 


27r 


Im {p*(r)q(r)} , 


(B.6) 
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where A := 2ir/k 0 is the wavelength in the background medium. 

Our equations are obtained by means of the “direct method”, i.e., from the representation 
formulae (essentially Green’s second identities) applied to boundary values of the fields in all 
material regions. Such a formula for a field p satisfying Eq.(B.2) in a region is 


L/ r i^wr p(r,) -^ SJr - r ' )q(r,) . 

[ dV [d' n g m (r - r') p(r') - g m (r - r ') q{r')\ = 
JdQ 


p(r) for r e Cl m , 
0 for r qL f2 m , 


(B.7) 


where the derivatives are taken along the normal exterior to Q m (as in Fig. 15), p and q in 
the integrand are the boundary values of the fields, and 


9 m (r) = 




47t r 


(B.8) 


is the Green function of the Helmholtz equation in the region h?. 

While the formula (B.7) holds in the entire space, integral equations relate boundary values 
of the fields. They can be expressed in terms of the single- and double-layer boundary operators, 
belonging to the standard set of operators: 



■=— d 2 r' g m {r-r')4>{r') 

(single-layer) , 

(B.9a) 

PQ Jr m 



(K m <P)(r) 

m 

(double-layer) , 

(B.9b) 

«,«>)(<•) 

■■=[ dV 0s "( r ~ r V) 

Jr dn(r) 

m 

(adjoint double-layer) , 

(B.9c) 


Pm dn(r) Jr m dn(r') 

(hyper-singular) 

(B.9d) 


(our definition does not include the factor “2” appearing in some references (e.g., in [10]). 

[we do not enter into a rigorous discussion of the conditions under which the representation 
formulae are valid, including the question of the functional spaces of the sources and fields - 
this is a complex subject; we only mention that most of the results in the following hold not 
only for smooth surfaces, but also for Lipschitz-continuous surfaces, hence, in particular, for 
polyhedral domain boundaries commonly assumed in practical applications 

In the notation of Fig. 15, the boundary values of the integrals in the representation formula 
(B.7) are now 


and 



^ r ' d 'n 9 m {r - r') p(r') 

= ( K m p)( r )±hp( T ) 

for 

(B.lOa) 

r d 2 r , g m (r — r') q(r') 

= (V m p)(r) 

for 

(B.lOb) 


m 
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Hence, for the observation point approaching the boundary from either side, the representation 
formula (B.7) yields the relation 


b= - K m P+ V m<l ■ (B- 11 ) 

By multiplying the above equality by Po/Pm? faking its normal derivative, and using K' m = 
(Po/Prn) d n v m and W m = ~ (Po/Pm) d n K m > we find another relation, 

\ ( l = w m P + K l m q . (B.12) 


Analogous relations are obtained for the exterior region i? 0 , where the representation formula 
includes also the incident fields (asymptotic field values on the “boundary” of i? 0 at infinity. 


The resulting set of four relations is then 

W 0 p+{±I + K , 0 )q = qW , (B.13a) 

(±I-K Q )p+V 0 q=p<>™) , (B.13b) 

-W m p+{\I-K' m )q = 0, (B.13c) 

(\I+K m )p-V m q = 0. (B.13d) 


By taking the difference of (B.13a) and (B.13c), and the difference of (B.13d) and (B.13b), 
one obtains the manifestly symmetric first-kind system of equations 


'^0 + ^m 

K + K' m 


p 


g( inc) 

. K 0 + K m 

-v 0 -v m _ 


_q_ 


— p( inc ) 


which can be proven to be strongly elliptic [11]. Similarly, by taking just equations (B.13b) 
and (B.13d), one finds the “simple equations” for the transmission problem [12, 13, 14, 15], 

_\I + K m -V r 

Further, taking the sum and the difference of the above equations, results in the set of 
equations 

'l~K 0 + K m V 0 -V, 

.~K 0 -K m V 0 + V r 

which, due to the appearance of the identity operator in one of the diagonal blocks, can be 
considered a second-kind system. 



P 


p { inc ) 


_q_ 


0 


(B.16) 



P 


pi inc ) 


_q_ 


0 


(B.15) 


B.2.3 Surface integral equations for p(r) and dp(r)/dn(r) for multiple regions and 
interfaces 

In this Section we present a brief derivation of the surface integral equation formulation for 
problems involving multiple material regions. The two unknown functions appearing in such 
integral equations are chosen to be the pressure p(r) and its normal derivative, both defined 
on defined on material interfaces. 


31 




















Although we have obtained such equation in both first- and second-kind formulations, we 
discuss here, for definiteness, the system of second-kind integral equations - a generalization 
of the system of equations (B.16). 




Figure 16: Examples of simplest multi-region geometries: nested regions (a) and adjacent 
regions with a triple-interface junction (b). 

In a multi-region transmission problem we consider regions fl m and oriented interfaces r i with 
normals pointing from their negative to the positive sides. We introduce “orientation factors” 
a im = ±1, defined such that 


O'A 


+ 1 if fi m is adjacent to the positive side of F i , 
< — 1 if fi m is adjacent to the negative side of F { , 
0 otherwise . 

V 


(B.17) 


In this notation, the relations (B.13b) and (B.13d) generalize to 

E(Fi jPj - a jm K mPj + a jm V m<lj) = S m0 Pi^ ( B ' 18 ) 

3 


for any interface r i belonging to the boundary and for any orientations of interfaces. 

By considering the other region adjacent to say i?^, we obtain the analogous relation 

S ikPk ~ a kn K nPk + a kn V n %) = 6 n0Pi^ > ( R19 ) 

k 

where the sum is taken over interfaces r k forming the boundary of fl n (Fig. 54). 
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Figure 17: An interface F i and the two adjacent regions with their boundaries. 
In the matrix form, the system (B.18) and (B.19) can be represented as 


\l-a- K cr. V 


Pi 

2 im m im m 


2 I ® in ^n ® in ^n 


y. 


+ £ 


®jm ^m ® im ^m 


Pi 


+ £ 

k^i 


® kn ^n ® kn ^n 


- 



/ nc) " 

°mOPi 


Pk 

= 

s » (inc) 

°n0 Pi 


%_ 



(B.20) 


By taking the sum of the first and second rows of these equations and the sum ot the rows 
multiplied by a im and <J in , and by rearranging the terms, we obtain, for each interface F ^ a 
system of equations of the form 


LA, 


- 


2 (inc) 

a i0Pi 

Pj 

= 

(inc) 

fioPi 




(B.21) 


where the blocks of operators are given by 


Aj - s zj 


I 

0 


-cr. K a- V 



+ E 

jm m jm m 

0 

0 

— a■ a ■ K a- a ■ V 



m 

im jm m im jm m 


(B.22) 


where the sum runs over all regions Q m shared by both the interfaces F i and F ■. In the case 
of a single material domain, the system (B.21) reduces to the sum-and-difference equations 
(B.16). 


B.2.4 Scaling of matrix blocks for high-contrast problems 

In high density-contrast problems equations (B.16) and their multi-region generalization (B.21) 
exhibit unfavorable conditioning properties: Since the single-layer operators V m (Eq. (B.9a)) 
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are proportional to the factors p m /p 0 1 (for m 7 ^ 0 ), the operators V m dominate V 0 
(wherever the latter appear) and all blocks involving single-layer operators are proportional 
to large density ratios. For instance, at low frequencies, the block of operators in the system 
(B.16) becomes 


'l~K 0 + K m 

O 

1 

3 

_1 


1 

-Vm' 

. ~ K o~ K m 

V) + Vn. 


. - 2K 0 

v m _ 


we used here the fact that, for an object much smaller than the wavelength, there is little 
difference between the operators K 0 and K m (provided the refraction indices of the media are 
similar, which we always assume). Hence, for a high density contrast, there is an imbalance 
between the operator blocks and the condition number of the system increases by a factor of 
order p m /p Q . 

A simple way of improving the conditioning of the system is to multiply all the matrix 
blocks involving the single-layer operators by a small factor £ of order Po/Pm? s °l ve the modified 
system, and finally multiply the q blocks of the solution by the same factor £. For instance, 
after solving the rescaled system (B.16), 


I-K Q + K m 

( (v„ - VJ ' 


p 


^(inc) 

-*0 ~ K m 

t(V 0 + VJ_ 


A. 


0 


(B.24) 


the actual solution is obtained as q = £ q. 

As we discuss in Section B.5.1, the above preconditioning procedure improves very signif¬ 
icantly the convergence of the iterative solution. More importantly, for low frequencies and 
larger contrasts (p m /po ^ 100 ), preconditioning is necessary in order to obtain the correct 
solution of the system, unless one is prepared to carry out a large number of iterations to reach 
relative residual values of order 10 -6 or so; and even in that case the quality of the solution is 
questionable. 
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B.2.5 Coupled volume and surface equations for partly inhomogeneous objects 

We present here a formulation of a system of integral equations describing an object consisting 
of both homogeneous and inhomogeneous regions; this type of a problem is illustrated with 
two examples: 

A homogeneous inclusion in an inhomogeneous object. We consider an inhomoge¬ 
neous region Q embedded in background-space region i? 0 and containing a homogeneous in¬ 
clusion fi m (Fig. 18). 



Figure 18: An inhomogeneous region Q with a homogeneous inclusion f2 m . Symbols and 
r± in the text refer to positive and negative sides of the interfaces, according to the directions 
of the normals. 


In the region M 3 \ f2 m the Euler equation satisfied by the pressure p can be written in an 
equivalent form 


(V 2 + kl) p(r) - kl (l - ^) p(r) - V • 
=: (V 2 + kl) p(r) + S p (r ) = 0 



Vp(r) 


(B.25) 


where the source term S p is a functional of p (dependent on p and Vp). 

It follows from Eq. (B.25) that the pressure field satisfies the representation formula 




' dgoCr-r') 
dn(r') 
d%(r ~ r') 
dn(r') 


P{r’) - 9 0 (r - r 1 ) 

P(r') ~ ffo( r ~ r> ) 


dp(r') 

dn{r r ) 

dp(r') 

dn(r') 


[ d 3 r'g 0 (r 
JQ 


' 5,(0 



for r ^ 12 , 
for r E 12 . 


(B.26) 
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The two surface integrals in the above equation are due to the outer and inner parts of the 
boundary of 17, dQ = T 0 U (Fig. 18). The derivatives are taken, respectively, along the 
normals exterior to the domains 17 and l7 m , as also shown in Fig. 18. The field and field 
derivative are the limiting values reached when approaching the boundaries from the interior 
of the domain 17. 

In order to obtain the L-S equation in 17 we also consider the scattered field, satisfying the 
representation formula 


d% (r - r') 
dn(r f ) 


P {sc) (r')-g 0 (r-r') 


d p( sc )(r') 
dn(r') 



for r G 17 0 , 
for r 0 17 0 . 


(B.27) 


The normal direction in Eq. (B.27) is the same as in (B.26), i.e., exterior to the domain 17 (as 
shown in Fig. 18) and the fields in the integral are taken on the “interior” side of the boundary 
of 17 0 bounded by a large sphere of radius i?, df2 0 R — r 0 U T R . The Sommerfeld radiation 
condition implies that the contribution to the integral from r R can be dropped. 

Similarly, the solution p(r) of the Helmholtz equation (B.2) satisfies the representation 
formula 


dV 


S<, m (r-r’) p (r ') _ 9m(r _ r ' ) 9p ( r ') 


dn(r') 


dn(r') 


0 for r fl m , 

— p(r) for r <E Q , 


(B.28) 


with the normal n m shown in Fig. 18 and with the Green’s function in the region J? m , 


9 m (r) = 


e ik m r 

4i r r 


(B.29) 


The representation formula (B.28) can be obtained in the same way as (B.26) in the absence 
of the source term S. Further, by using the boundary conditions (B.l) to relate the values and 
gradients of p taken on the two sides of the boundary T m , we can rewrite the representation 
formula (B.28) as 



-Q 

dn{r f ) 


p( r ') — 9m( r — r ') 


Pm dpjr 1 ) 
p{r f ) dn(r') 


0 for r 0 I7 m , 

-p(r) for r G I7 m , 


(B.30) 


where p(r'), p(r'), and d n p(r f ) are the limiting values of these functions for 17 3 r' —T m . 

By considering r G 17, taking the difference of Eqs. (B.26) and (B.27), and using the 
relation p — p ( sc ) = p( mc ) we obtain 


L 

d L rr 


d V 


+ / d 2 / 


' 9g 0 (r-r') 

dn(r') 

' dgoir-r 1 ) 

dn(r’) 


P(r') ~ 9 0 (r - r') 


z dp{r' 


dn(r') 


[ d 3 r' g Q {r -r')S {r') 

JQ 


P(r') ~ g 0 (r - r') dp(yT ^ 


dn(r') 


(B.31) 




i 2 / 

d r 


^ (r-o (sc) _ dp^y 

dn{r’) 0 dn(r') 


= — p(r ) for r E Q . 
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By using again in the last integral p ( sc ) = p — p( inc ^ we find 

- L dV [ %' } p<r ' ) ~ 9 ° (r ~ r ° ~ X dVa « (r ■ r ' )s '’ <r ' ) 


+ L~L ' dV ^ - -' ) ap<r,) 


L 


+ / dr 1 


o J o 
2 / 


dn(r') 


dn(r') 


(B.32) 


dgo( r ~ r/ ) ,(inc) (r / ) _ (r _ ^ 9 ^ (lnC) (^) 
dn(r') 0 dn(r’) 


— — p(r ) for r G fi 


But the last integral in this expression is, by the Green’s second identity, simply — p( mc )(r). In 
the previous difference of integrals the terms proportional to p cancel, because p is continuous 
across T 0 - However, if p is discontinuous across the boundary, dp/dn is discontinuous as well. 
The boundary condition (B.3) yields then 


/, 

/, 


d v 


dg 0 (r-r r 


7 ^- P(r ') -g 0 (r- r’) d 3 r’g 0 (r - r') S p (r') 


dn(r 


d r'g 0 (r 


Po 


dp{r') 
p(r f ) J dn{r f ) 


_p( mc )( r ) = —p[r^ for r G 17 . 


(B.33) 


The above expression, with the source density S p given by Eq. (B.25), constitutes the volumet¬ 
ric equation in the desired system. The surface equation is simply the representation formula 
(B.30) with the point r approaching the positive side of T m . Thus, finally, the set of coupled 
volume and surface equations is 


p( mc )(r) = p(r ) 


ft(r') 


+ f^ r ' g 0 ( r - r- 1 ) l k% 

-/,/ 


p(r') + V r , 


v{r ') _ Q ( r _ r '\ dKO 

dn(r') 0 dn{r') 


and 


d 2 r 


r+ 

A m 


dn(r f ) 


P(r') ~ 9m( r ~ r> ) 


/\ Pm dp(r') 


p(r') dn(r') 




for r G h? 


= 0 for r G TA . 


(B.34a) 

(B.34b) 


In Eq. (B.34a) we marked explicitly that it holds not only inside the domain h? but, by 
continuity, also on its boundary. 


Adjacent homogeneous and inhomogeneous regions. As another example we consider 
adjacent inhomogeneous and homogeneous regions Q and i2 m , both embedded in background- 
space region f2 0 and containing (Fig. 19). 
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Figure 19: A homogeneous region ft m adjacent to an inhomogeneous region ft. T m and T m0 
denote two parts of the boundary of ft m . 


The representation formula (B.26) remains in this case unchanged, provided T+ is the 
interface of ft and ft rn . On the other hand, the analogue of the formula (B.27) for the scattered 
field contains an additional contribution from the interface between the regions ft 0 and i? m , 


d V 


T+urj 


-r' 

on(r') 


dp( sc \r') 

dn(r') 



for r £ ft 0 , 
for r ^ ft 0 , 


(B.35) 


where we dropped the vanishing contribution from r R . By taking again the difference of 
Eqs. (B.26) and (B.35) for a point r G ft and by using p ( sc ) — p — we obtain 


[ d2r ' dg ^ V ( ^ ^(0 ~ g ^ r ~ r ') ~ [ ^ r ' 9 ^ r ~ r ')S (r') 

Jr+ur+ 0 L dn(r>) dn{r')\ J Q p 

\L-L 




+ 4 dV 


v C™)(r') -a (r- r') d W 
dn{r') 0 dn{r') 


dn{r r ) 

= —p(r) for r G ft 


(B.36) 


which differs from Eq. (B.32) only by the additional contribution from the interface 0 . The 
remaining steps in the derivation are as before, and the final set of coupled volume and surface 
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equations becomes 


p( mc )(r) = p(r) 

+ f^r' g 0 ( r - r> ) { k o 


1 - 


K (r') 


-/ dV 9 o (r-,') (l-^y) 


Po ^ 

dn(r') 


p( r ') + V_, 




dV 


/ r+i jr+ 

J ™ Ui m 0 


a go( r ~ r/ ) p(r /) _ q ( r _ r /) Ml' 

dn(r') 0 dn(r') 


and 

L 


r+\ j r+ 

1 m^-L rn Q 




<9n(r') 


p(r') dn(r')_ 


for r G i? (B.37a) 

= 0 for r € T+ U T+ 0 , 

(B.37b) 


where the surface equation (B.37b) is simply the representation formula for the region Q m 
involving the entire boundary of that region. 


Discretization. For concreteness we consider here the geometry of Fig. 18 and write the 
resulting system of equations (B.34) as 

p{r) - (U Pq) (r) 

- (£ 0 Pr+)W + ( v o d nPr+)( r ) =P {mc \r ) for r € f2 , (B.38a) 

((P + KJp rj; ){r)-(v m ^d n p ri )(r) = 0 href.. (B,38b) 

Here in Eq. (B.37a) the operator U is a generalization of the volumetric “Newton potential”, 

{UPn)( r ) ■= j^r’ g Q (r - r') S p {r') , (B.39) 

expressed in terms of the source density (as defined by Eq. (B.25)); V 0 is the single-layer 
potential 

(r 0 4>r)( r ) ■= J^d 2 r r g 0 (r - r') cf>(r r ) , (B.40) 

due to the source (j) on an interface F, and JC 0 is the analogous double-layer potential, both 
defined at all points of the domain Q. The operators K m and V m in Eq. (B.37b) are the usual 
double- and single-layer boundary operators, the latter not including the density factor present 
in the definition (B.9a). For simplicity, we assumed here material properties are smooth on 
the boundary F 0 , which allows us to omit the contribution of this surface to the volumetric 
equation (the third line of Eq. B.34a). The presence of absence of such contributions (related 
to discontinuities of the material parameters [9, 1, 2]) is a problem belonging to the purely 
volumetric formulation, unrelated to the coupling of volume and surface equations. 

With F+ denoted simply by F, the unknown fields appearing in Eqs. (B.38) are the volume 
pressure p^, surface pressure p r , and its normal derivative d n p r . For simplicity of notation, 
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we assume now that Pq is expanded in constant trial functions <p t supported on tetrahedra t of 
17, while the p r and p r are represented in terms of constant functions supported on facets 
/ of r. Accordingly, we denote the vector containing the unknown expansion coefficients as 
consisting of three blocks, [p t d n Pf ] T . 

In devising the testing procedure, we recall that Eq. (B.38a) holds not only in the region 
17, but also on its boundary T. Therefore, that equation should be projected on both volume 
and surface testing functions, for which we again assume the functions p t and </>j. Then, when 
projected on the latter functions, the potential operators /C 0 and V 0 become the boundary 
operators K 0 and E 0 , and in this way one recovers the full set of surface equations. The 
resulting discretized equations take then the schematic form 


-(£»)« ( v »).f 


Pt 


(inc) 

Pt 

('-«) * o)« (V'o)ff 


Pf 

= 

(inc) 

Pf 

0 (\l + K m ) ff -{p/p m V m ) ff 


d nP{ 


0 


(B-41) 


with a square coefficient matrix. 

In the absence of the homogeneous region l7 m the system (B.41) reduces trivially to the 
purely volumetric L-S equation. It if also easy to show that, when region 17 becomes the 
background space, Eqs. (B.41) reduce to the usual surface equations (B.15) for the transmission 
problem, which can be then recast into the sum-and-difference form (B.16). In general, the 
system (B.41) can be similarly rearranged by taking the sum and the difference of the second 
and third equation. It is also evident that other than constant basis functions can be used in 
the discretization. 
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B. 3 Implement at ion 

In this Section we describe some freatures of the current implementation of the solver, in 
particular discretization of the integral equations for multiple material regions and compression 
of the resulting coefficient matrix. 

B.3.1 Aspects of the discretization 

Choice of trial and testing functions. The approach often used in discretizing trans¬ 
mission problems (involving single- and double-layer operators, and possibly hyper-singular 
operators) is to discretize the pressure and velocity fields p and q with piecewise constant and 
piecewise linear functions, respectively. This choice is motivated by the mapping properties of 
the operators [11] appearing in Eqs. (B.14) and (B.16): the operators K are of order 0, while 
V’s and W’ s have have orders —1 and +1. 

From the point of view of the smoothing and differentiating properties of the operators, a 
possible drawback of a linear interpolation for q is that it enforces continuity even at sharp 
edges of the surface, where q is expected to be discontinuous. However, in anatomical models 
we are working with the surfaces tend to be relatively smooth, linear interpolation for both p 
and q appears to be advantageous. In fact, in our code we have implemented linear-functions 
discretization for p and both constant- and linear-functions discretizations for q and our ex¬ 
perience has shown that, overall, linear node-based functions are preferable. With this 
choice, sharp edges and junctions of the surfaces can be handled by simply relaxing continuity 
constraints on the velocity field, i.e., by introducing independent variables on edges of separate 
surfaces. 

As to the choice of testing functions, most engineering applications use collocation, in order 
to minimize the cost of the stiffness matrix construction (as discussed, e.g., in [16]). In our 
implementation we opted for a Galerkin method with the same spaces of testing and trial 
functions; its higher computational cost is offset by matrix compression which requires only 
computation of a highly sparse matrix representing couplings of nearby elements. 

Treatment of junctions of surfaces. We first describe our treatment of junctions of tree 
interfaces in the context of the previously discussed system of a single body partitioned into 
two adjacent regions, shown in Fig. 16(b). 

In this case the non-vanishing orientation parameters can be written as 


(7 u ~ ~ 1 5 a io = +1 
a 22 = — 1 5 a 20 = +1 
(J32 = —1 , cr 31 = +1 


(B.42a) 

(B.42b) 

(B.42c) 
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The final matrix of operators is, therefore, 


I-K. + K, 

Vo-Vx 

o 

1 

V 


Vx 

-k 0 -k x 

v+v 

o 

1 

V 

Kx 

~Vx 

o 

1 

V 

i-k 0 + k 2 

V 0 -V 2 

k 2 

-v 2 

1 

o 

V 

-k 0 -k 2 

v 0 + v 2 

1 

to 



-v 


-v 

i-k 1 + k 2 

v-v 

Kx 

~Vx 

-k 2 


-Kx-K 2 

t>i + t> 2 


and the full system of equations takes the form 


Pi 


(inc) 

Pi 

r 


(inc) 

Pi 

P2 


(inc) 

P 2 

$2 


(inc) 

P 2 

Ps 


0 

% 


0 


(B.43) 


(B.44) 


According to the continuity criteria imposed on the solutions, we to discretize the system 
(B.44) in terms of the total of seven sets of “composite basis functions” (CBFs), by expanding 
the unknown fields in trial basis functions as 

p(r) = ^2p vl $ vl (r) + ^Pv 2 $A r ) + XX3 ^v 3 ( r ) + XXj T vj( r ) 

vl v2 v3 vJ 

and 

Q( r ) = ®V1 ®Vl( r ) + 5^ ^V2 ^V2( r ) + ^V3 ^V3( r ) * 

VI V2 V3 

All the basis functions are linear functions associated with vertices (nodes) of the surface 
mesh: <P vi are supported on sets of facets sharing interior vertices of an interface F ^ F vi are 
supported on sets of facets sharing all vertices (including boundary ones) of the interface, and 
T vJ are associated with vertices of a junction J, i.e., with common vertices of several (here 
three) interfaces. 


(B.45a) 

(B.45b) 
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Assuming the same as above set of testing functions, we arrive at the discretized system 
described by the square matrix which can be represented as 


^vlvl 

M vlVl 

^vlv2 

^vl V2 

^vl v3 

^vl V3 

^vlvJ 

Mvi vl 

-^vivi 

Mvi v2 

^V1V2 

^V1 v3 

-^Vl V3 

-^VlvJ 

-^v2vl 

-^v2Vl 

Vv2 v2 

-^v2V2 

-^v2v3 

-^v2V3 

-^v2vj 

M V2vl 

M V2V1 

M V2v2 

M W2W2 

M V2v3 

M \2\3 

M V2vJ 

-^v3vl 

-^v3Vl 

-^v3v2 

M v3V2 

-^v3v3 

-^v3V3 

-^v3vJ 

^V3vl 

-^V3V1 

-^V3v2 

^V3V2 

-^V3v3 

^V3V3 

^V3vJ 

^vJvl 

M vJVl 

-^vJv2 

M vJV2 

-^vJv3 

M vJV3 

^VJVJ 


Here the 2x2 blocks involving the “ordinary” basis functions and (associated with single 
interfaces) can be easily constructed from the operator matrix representation (B.43). For 
instance, 


and 


^v2v2 

-^v2V2 


\l-K a + K 2 )^ 2 

(£« - 

Mv2v2 

^V2V2_ 


. ^0~ -^2)v2v2 
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-^v2V3 


(-^2)v2v3 

( ^V)v2 V3 

_-^V2v3 

-^V2V3_ 


.( — -^2)v2v3 

(^2)v2V3 . 


(B.47) 


(B.48) 


etc. 

The blocks in the last row and the last column, involving the junction basis functions, have 
a more complex structure. To exhibit them explicitly, we have to represent the basis functions 
T vJ (still schematically) as 

3"vj = r VJ i + r VJ2 + r VJ3 , (B.49) 


where T vJi is supported on facets belonging to the interface r i . 

Suppose then we want to find the formula for the block M vJ vJ . We go back to the expression 
(B.22) and take its (1,1) block involving matrix elements between basis functions representing 
the pressure, i.e., 5^ I + J^ m (— a jm^m)^ an d obtain 


E E( r v« • (*« 1 - E <v K J r vj j) 

i=1 j=1 m 

= E(r VJi . ir *n) -EEEv . 

i=l i=l j=1 m 


(B.50) 


This expression is still schematic, as it represents an entire matrix block. To be more specific, 
we may consider the matrix element between the basis functions T vJia and where a 

and /3 label vertices on the junction (Fig. 20). 
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Figure 20: One of contributions to a matrix element between two junction-type basis functions 
associated with vertices a and (3. The arrows indicate orientations of interfaces and the dashed 
line shows the coupling of the facets f i cT i and fj C r ■ through the operator K x in the 
region f2 1 . 

The contribution indicated in this Figure is the second term in the last expression of Eq. (B.50) 
with i = 1, j = 3, and m — 1, i.e., a coupling of facets /q and / 3 on interfaces r 1 and F 3 
through the region £2 1 . 

In the actual implementation in the code, the junction-junction matrix elements such as 
given by Eq. (B.50) are evaluated together with all other matrix elements between ordinary 
basis functions. In fact, the considered facets f 1 and / 3 are also members of ordinary basis 
functions ^ vl and ^ v3 supported on the interfaces r i and F 3 , and the contributions to various 
basis functions are being evaluated simultaneously. 

More specifically, our matrix fill algorithm involves an outer loop through pairs of interfaces 
r i and r '• and an inner loop through pairs of facets f i C and fj C r In the latter loop 
all matrix elements between “primary” basis functions supported on f i and fj are computed 
(there are 3-3 = 9 such matrix elements, since each triangular facet supports thee different 
linear basis functions associated with three vertices of the triangle); then, contributions of facet- 
facet matrix elements are being added to all relevant matrix elements between the “composite” 
basis functions, including junction-type functions. In other words, the formula (B.50) is not 
implemented literally as the sum indicated there, but rather by evaluating that expression 
“inside-out”: we first compute the matrix elements (T vJ ^ a , K m Y v jj^ for all regions 
and for all pairs of “primary” basis functions labeled by interfaces r i and facets f i and 
fj, and vertices a and /?, and then add their contributions to the matrix elements in all blocks 
in (B.46). 

B.3.2 Matrix compression and fast matrix-vector multiplication 

In the approach we describe here, we utilize the FFT-based AIM compression technique [5] 
initially developed in the context of electromagnetics and adapted to acoustics [1]. Such a 
matrix compression was developed in order to enable solving large linear sets of equations 
with dense matrices utilizing storage and execution times characteristic of problems involving 
sparse linear systems. (The physical idea behind the compression methods is that interactions 
at large distances require less resolution than interactions at small distances. As the result, 
the computational complexity and memory requirements of the compression methods scale 
approximately linearly with the number of unknowns N.) 

In order to explain the structure of this algorithm on a simple example, let us suppose 
we are evaluating the far-held contribution p f to the pressure, generated by the double-layer 
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operator K acting on the input pressure field. We assume the input pressure has the expansion 

p{r) = ^ j p a T a (r) = y ^p a ^ j T i a (r) . (B.51) 

OL Oi i 

Here as are the unknown numbers and the index i labels interfaces r\ supporting the composite 
basis functions (CBFs) T^. Hence, in general, the sum 

rjr) :='£rZ(r) (B.52) 

i 

constitutes a junction-type basis function supported on several interfaces; a particular case are 
“ordinary” basis functions associated with single geometry components. 

We are the evaluating projections of Kp on the same set of basis functions, i.e., 

p' a :=(T a ,Kp) . (B.53) 


For the purpose of evaluating the far field, the basis functions T a are approximated by distri¬ 
butions of sources on nodes of a Cartesian grid, 


T a( r ) A au S3 ( r ~ u ) 


(B.54) 


With the operator kernel K(r,r') = dg(r — r')/dn(r '), the result of the operator acting 
on the field is 


(■ Kp)(r ) = J2 p 0 f d3 


, dg(r - r' 


dn{r 


7 y Lr ^ , ) = Y,P^ / dV 5( r - r ')^)V). 

(B.55) 


where the Green function g is taken with the parameters of the region through which the points 
r and r' communicate. 

Suppose first the equations involve only a single operator 1C of the form such as one of the 
blocks of (B.22). In this case, the far-held computation procedure is as follows: 

We assume that the source x is represented in terms of the trial basis functions T x 

x \ r ) = y }2 x i3 T 0 l ( r ) = T f J ( r ) • (B.56) 

P P 3 

Here as are the unknown numbers and the index i labels interfaces F i supporting the composite 
basis functions (CBFs) X T^. We then consider the held y = /C x created by an operator /C acting 
on the sources, projected on the testing basis functions T^ z , 

y a :={TZ,Kx) . (B.57) 

The operator is assumed to be of the form 

JC = D y gD x , (B.58) 
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where g is the Green function of the region in which the operator acts, and D y and D x are 
differential operators (e.g., normal derivatives) acting at the field and sources points, hence 
the field projections can be, after integrating by parts, represented as 

Va = J 2 [ d2r i / d 2 r 2 D y T a k ( r i ) 9 ( r ! - r 2 ) D x T^ l (r 2 ) x l 0 . (B.59) 


At large distance between the field (observation) and the source points the (derivatives of) 
the basis functions appearing in Eq. (B.59) can be approximated, for each interface F ^ by 
distributions of sources on nodes of a Cartesian grid, 


D T y 

± ol 


k i 


(v) X J2 A ( D y r ^')^ 3 (r-u) , D x r* H (r) *'£ A ( D * r °‘)v S3 ( r -v) ■ (B.60) 


The basis function is then approximated as 


0/(0 - J2 A f,q 63 ( r ’ 


(B.61) 


One of possible expressions for the far-field contribution to the coefficient matrix results from 
substituting the expansions of auxiliary basis function given in terms of monopole sources is 

A t.h = E ■ (B ' 62) 

Ql 


In the case of the general integral-equation problem involving a matrix of operators and 
the corresponding vectors of solution and r.h.s. functions, as exemplified by Eq. (B.44), the 
far-field contribution to a matrix block acting on a source vector x has the form 


sju - v) °U A ( D * 4 

m i u,v j l (3 

E EE d ( D y T a kt )u a im J2Sm( U ~ V ) E a jm A ( D x T ^' J )v 4 

m iu v jig 


(B.63) 

= E E A ( D y E s™(“ - ®) V 

m iu v 

s -V- / 

Y m 
1 u 

m iu 

The formula (B.63) is implemented as an outer loop through regions i? m , containing nested 
loops through the indices (j, /?), v, and (i,u). The first inner (j,/3)-loop evaluates the “Carte¬ 
sian vector” X™, i.e., the equivalent-source representation, in the region i? m , of the input 
MoM vector of coefficients Xp. The v-loop is implemented in terms of FFTs: the Cartesian 
source vector X™ is transformed to the Fourier space, multiplied by the Fourier transform of 
the Green function g m , and transformed back to form the Cartesian field vector Y™. In the 
loop that vector is converted to the output MoM representation of the vector y. All the 
above operations are performed for all regions f2 m (the outer loop) and the contributions to 
the far output field are being added. 
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B.3.3 Parallelization 


In view of the relatively large scale of the problems we are solving - a reasonably realistic 
geometry requires at least several hundred thousand unknowns - not only matrix compression, 
but also parallelization of the solver becomes a necessity. 

In our previously developed volumetric solver [2, 1] we have implemented a distributed- 
memory message-passing (MPI) parallelization in conjunction with the FFT-based matrix 
compression. It involved partition of the object geometry into multiple slice-type segments 
and assigning each of the slices (together with its two or more nearest neighbors) to individual 
processors. Each processor was using its “main” slice and its neighbors to construct a geometry 
of a size sufficient to evaluate all the near-field couplings of the main slice elements with 
the remainder of the object. The far-held couplings were handled by means of the matrix 
compression and fast matrix-vector multiplication scheme, which also required each processor 
to store and handle only quantities dependent on its local geometry and only the corresponding 
segment (slice) of the global Cartesian grid. 

One of the main reasons that distributed-memory approach has been effective was that the 
geometry was volumetric and relatively uniformly discretized . These circumstances contributed 
to a reasonable work load balancing, while the regular structure of the geometry partition 
allowed a fairly simple communication scheme in exchanging data between the processors and 
their memories. 

In contrast, our present solver has to be able to handle surface models with fine geometry 
details and widely varying discretization scales. In fact, in our models the triangle sizes span 
two order of magnitude, from about 5 mm in smooth geometry regions to 0.05 mm = 50 /am in 
the middle and inner ear details. As a result, balancing the computational load by means of 
geometry partition becomes increasingly difficult. Even if the numbers of geometry elements 
(triangles, vertices) or unknowns were approximately equalized between the processors, the 
numbers of their near-field couplings (and thus matrix elements) would still vary in a wide 
range, typically from few tens to several thousand. 

In this situation we opted, at least as an interim solution, for implementing a shared- 
memory parallelization, especially in view of the presently available computational hardware. 
Current desktop or even laptop systems provide relatively large amounts of RAM storage (tens 
of gigabytes) accessed by several processors. Higher numbers of processors, along with larger 
storage, are available on supercomputer architectures, typically allowing both message-passing 
and shared memory access. 

Our experience has shown that the shared memory (OpenMP) parallelization of the present 
solver is entirely feasible and efficient. In the following we describe parallelization of the main 
operations in the solver, in the order of decreasing computational cost: 

1. In our nonuniform discretization problems the most computationally intensive solution 
stage is the “far-field subtraction” in stiffness matrix construction. This operation (an 
essential part of the AIM algorithm) handles the precomputed near-field matrix and 
subtracts from them the corresponding far-field elements. Although the near-field is 
globally sparse, its parts associated with finely discretized geometry regions may be 
quite “dense”, i.e., matrix rows may contain several thousand nonzero elements. The 
main computational cost is, therefore, evaluation of the far field (expressed in terms 
of the equivalent-source coefficients and the Green function tabulated at nodes of the 
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Cartesian grid) for each of the near-field matrix elements. At the same time, however, 
this operation is performed independently for each matrix element, hence the resulting 
loop can be immediately distributed among the threads. As a result, a nearly perfect 
parallel speedup can be achieved. 

2. The next most computationally expensive task is the evaluation of the near-field stiffness 
matrix itself, chronologically preceding the previously described subtraction step. It 
consists of two stages: 

(a) Computation of the sparsity pattern of the matrix, based on the near-field range. 
In this stage we generate a special structure (not yet the stiffness matrix itself) 
storing information on the elements to be included in the matrix. Depending on 
the type of basis functions (e.g., those associated with facets or with vertices), the 
sparsity pattern stores, for each facet or vertex, the list of near facets or vertices. 
The relevant geometry elements (facets or vertices) are here first distributed among 
geometry buckets, each assigned to one and only one bucket. Then a double loop 
through pairs of nonempty and sufficiently near buckets is executed and, within 
those loops, loops through pairs of geometry elements in the buckets. In the latter 
loops the distance the elements is checked and, if it does not exceed the near-field 
range, the appropriate element is stored in the sparsity pattern. 

Since the geometry elements are exclusively assigned to buckets, the outermost loop 
through the buckets can be split into threads without causing storage conflicts: each 
row of the sparsity pattern is then handled by only one thread. This scheme results, 
therefore, in a speedup nearly equal the number of threads. 

(b) Filling the determined structure of the sparse matrix. This task is more difficult to 
parallelize, 1 since the outer loops run through facets and facet-facet contributions to 
vertex-vertex matrix elements are being gradually added (this scheme is most effi¬ 
cient from the point of view of the total number of operations). Now, since different 
pairs of facets contribute to the same vertex-vertex matrix element, different threads 
may simultaneously try to write to the same storage location - although probability 
of such an event decreases with the problem size. In order to prevent such storage 
conflicts, we currently use the OpenMP locking function around the critical opera¬ 
tion of adding/storing contributions to the output matrix elements. This solution, 
however, is not fully satisfactory, as impairs the effectiveness of parallelization. A 
better scheme (in the process of implementation) is similar to that used in con¬ 
structing the sparsity pattern. The outer loops are taken through another set of 
buckets, now defined as containing sets of facets, but based on exclusive partition 
of vertices. This arrangement results in some redundancy: some facet-facet matrix 
elements are being computed several times, but only for facets near the boundaries 
of buckets. However, because of the exclusive assignment of vertices, it guarantees 
the absence of simultaneous store operations by different threads. 

Finally, the matrix-vector multiplication also involves two main operations: 

1 For basis functions associated with vertices; for basis functions supported on facets parallelization is trivial. 
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(a) Computation of the near-field component of the matrix-vector product. Paralleliza¬ 
tion of this task depends on whether the near-field matrix (stored in the sparse 
row form) is symmetric or not. In the latter case (which arises in most types of 
equations we are solving) parallelization is straightforward: the outer loop can be 
taken to run through the matrix rows, i.e., the storage indices of the output vector. 
Therefore, if this loop is split among the threads, each output address is handled 
by one thread only. 

There is no such simple solution in the case when the matrix is symmetric and only 
its independent (say, upper-triangular) part is stored. It seems that in this situation 
the optimal algorithm is to compute and store an additional list of matrix element 
numbers, organized by columns rather than by rows. This scheme allows computing 
the entire matrix-vector product in a loop through the output vector indices and, 
at the same time, increases the matrix size only by 1/3, rather than doubling it, as 
would be the case if the symmetric matrix were trivially stored as a nonsymmetric 
one. 

(b) Computation of the far-held contribution to the matrix-vector product. In this op¬ 
eration the main computational cost is the three-dimensional fast Fourier transform 
(FFT) performed on the sources or fields defined on the Cartesian grid covering the 
object. That FFT is implemented as a set of iterated one-dimensional transform 
taken in each of the three directions. The steps in the loops in which that transform 
is executed can be simply split among the threads without causing any conflicts. 


B.4 Some qualitative features of energy flow through high-contrast inter¬ 
faces 


Before discussing examples of computation, we give a short qualitative discussion of the ex¬ 
pected behavior of the energy flow in our problems. 

We first give a very elementary account of the acoustic reflection and transmission proper¬ 
ties in wave propagation through a half-space and through finite-thickness layers. We consider 
a one-dimensional problem of a layer of thickness a, filled with a material of density p and a 
refraction coefficient n, surrounded by a background medium of density p 0 and a unit refrac¬ 
tion coefficient. The layer is subject to a plane wave propagating in the direction normal to 
the its boundaries with the wave number k (in the background medium). The reflection and 
transmission coefficients (for the pressure field) are then 


and 




(1 — ( 2 ) sin (nka) 

(1 + £ 2 ) sin(nfca) + 2i£ cos (nka) 


(B.64a) 


T„ 


_2iC_ 

(1 + ( 2 ) sin (nka) + 2i£ cos (nka) 


e~ ika , 


(B.64b) 


where 

C := — (B.65) 

P 

is the acoustic impedance of the considered material. 
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For a “thin” layer, in the sense £ <C nka <C 1, the above expressions can be approximated 


by 


TZ a = l + 0((nka) 2 X 2 X/(nka)) , T a = 2i -j- [l + 0((nka) 2 , £ 2 , Q/(nka))] , (B.66) 


nka 

while for a large layer thickness (more precisely, for Im nka ^ 1, allowing a complex n), 
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(B.67) 


For comparison, for a semi-infinite half-space filled with the same considered material, the 
reflection and transmission coefficients are 


n = 


i-C 

i + c ’ 



(B.68) 


of the pressure wave entering the medium, is not small. 

The relative energy flux densities (Eqs. 
half-space problems are, respectively, 


K = 1 - l K J 2 = \rf = 


and 

F = 1 - \U\ 2 = Re £ |T| 2 = . (B.69b) 

Hence, the energy flux F a through the “thick” layer is (at least for real £) the square of the 
flux F entering a half-space medium - due to the fact that the wave is being reflected on two 
high-contrast interfaces. 


(B.4) or (B.6)) for the finite-thickness-layer and 


2C 


nka 


4C 


(1 + o 2 


(“thin” layer) , 


(“thick” layer) 


(B.69a) 




B.5 Examples of energy flow computations 

In this Section we discuss applications of the solver in modeling acoustic energy flow through 
the human head and its transfer to the cochlea. 

The most complete geometry model used consists of: 

(1) the outer surface of the skin surrounding the skull and containing 

(2) the outer ear represented by its exterior surface, 

(3) the surface of the auditory canal, 

(4) the tympanic membrane modeled as a finite-thickness surface; 

(5) the middle ear, consisting of the system of ossicles and supporting structures; 

(6) the skull, described by external surfaces of the bones constituting the skull and including 
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(7) a set of surfaces representing the inner ear cavity (boundaries of the cochlea, the vestibule, 
and the semi-circular canals). 

In order to assess the relative importance of various mechanisms of the energy transfer, 
we also discuss results of computations with modified or simplified models. In all geometry 
models we assume typical approximate values of the relative densities and the refractive indices, 
pj p 0 = 1000 and n = 0.2 for the soft tissues, and p/p 0 = 2000 and n = 0.4 for the bone. The 
model is being subject to a pressure plane wave of unit amplitude, incident from the left side 
of the head. The wavelengths A given in the following are always the wavelengths in the air. 

We first briefly characterize the models, and summarize the obtained results: 

(A) : A model of a homogeneous human head, filled with soft tissue. The computations in this 

case illustrate the significance of the scaling of the operator matrix and vector blocks 
in solving high-contrast transmission problems (Section B.2.4). We find that, in the 
absence of rescaling, although the residual in the iterative solution decays approximately 
exponentially, the slope is much smaller (in our examples 3 to 7 times, depending on 
the frequency) than the slope in the solution of rescaled equations. Furthermore, the 
solution of the non-rescaled equations converges to the correct solution only when very 
small residual values (10 -5 down to 10 -7 ) are reached, 2 and even then the resulting field 
distributions exhibit speckle-type fluctuations. 

(B) : A model of an isolated skull, surrounded by air. The computations, involving a rather 

complex skull structure, show the role of bone thickness in controlling the amount of 
the energy transmitted through two high-contrast interfaces, air-to-bone and bone-to- 
air. We also compare, for relatively low frequencies (wavelength A = lm) solutions 
of the transmission problem to the solutions of the hard-surface (Neumann) and the 
soft-surface (Dirichlet) problems; these problems serve as a verification of the soundness 
of our solution procedures. We find, in agreement with intuitive expectations, that 
the pressure in the Neumann problem tends to concentrate inside the skull cavity - in 
analogy to the expected behavior of the temperature distribution in the equivalent heat- 
conduction problem. On the other hand, the distribution of the velocity field in the 
Dirichlet problem is concentrated on protruding elements of the surface, in analogy to 
the charge distribution in the equivalent electrostatic problem. 

(C) : A model of a skull surrounded and filled with homogeneous soft tissue, with or without the 

outer ear canal. Here we can establish some properties of energy flow in a more realistic 
model with a rather intricate skull geometry, which includes the middle- and inner-ear 
cavity, and in particular the cavity housing the cochlea. We stress, however, that in our 
model we excluded the impedance-matching mechanism of the middle ear; hence, the 
energy entering the outer ear canal can only propagate through the surrounding tissues, 
and not through the usual air-conduction pathway. 

As a quantitative measure of the effect of the acoustic waves arriving at the inner ear 
through the head tissues, we compute the amount of energy flowing through the cochlea. 
More specifically, we evaluate the quantity <P S defined by Eq. (5.2) below as the average 
absolute value of the energy flux density through the surface of the cochlear cavity, 

2 Our computations are carried out in single precision. 
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relative to the energy flux density of the wave incident on the head. According to 
the estimate (B.69b) of the energy transmission through a high-contrast interface, the 
computed relative flux density is expected to be of the order of the ratio of the densities 
of the air and the tissues, i.e., ^ 10 -3 . 

The main results of the computations are that 

- The relative energy flux density <3? s is, as expected, of the order 10 -3 for large 
wavelengths (A larger than the head size) and, generally, decreases for smaller wave¬ 
lengths. 

- The presence of the outer ear canals increases the relative flux density by transport¬ 
ing the energy deeper inside the head. This enhancement is by about the factor of 
two at larger wavelengths (A > 30 cm) and the factor of three for 10 cm < A < 30 cm. 
For wavelengths A < 30 cm it exhibits a resonant behavior which, according to the 
computations with the this the next model (D) can be attributed to the outer ear 
canal. 

(D): A model of homogeneous soft-tissue head, without the outer ear canals and without 
the skull, but including the surface of the middle- and inner-ear cavity surfaces. We 
carried out computations with this model in order to assess the significance of the sound 
transmission through the bone vs. that through the soft tissues. As before, we compute 
the average relative flux density through the surface of the cochlear cavity. We find 
that for wavelengths A > 10 cm removing the skull from the head filled otherwise with 
soft tissue increases the flux by 10 % to 30 %. The results for smaller wavelengths indicate 
a resonant-type behavior similar to that seen in the absence of the skull and therefore 
attributable to the outer ear canal - the feature common to the two models. 

B.5.1 Exterior head surface (A): rescaling in solution of high-contrast transmis¬ 
sion problems 

As an illustration of the conditioning problems encountered in high-contrast problems (Sec¬ 
tion B.2.3), we discuss computations for a model of a homogeneous human head, assumed to 
be filled with the soft tissue. The model discretization gives rise to N — 58, 062 unknowns. 

In the first set of computations we assume the wavelength A = 1 m, corresponding to the 
frequency about 300 Hz. The resulting convergence histories for the original and rescaled sets 
of equations (Eqs. (B.16) and (B.24)) are shown in Fig. 21. 
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Figure 21: Convergence of the iterative solution to the transmission equations with the outer 
head surface, for A = 100 cm and the tissue density p m /p Q — 1000: no rescaling (£ = 1), and 
rescaling with the factor £ = 1/1000. 


Fig. 22 shows the energy flux distributions computed from solutions to the original and 
rescaled equations. Both solutions were obtained by iterating until the relative residual 10 -7 
was reached. At this very small residual value the energy fluxes are similar, but the distribution 
for the original equations exhibits distinct speckle-type irregularities. Further, the solution to 
the original equations approaches that for the rescaled equations only at this very low residual 
level, while the rescaled equations provide an accurate solution already for the relative residual 
of the order 10 -3 . 
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Figure 22: Distribution of the energy flux density in a model of a human head filled with a 
homogeneous tissue with p Yn j p 0 = 1000 and n m = 0.2, and for A = 100 cm, obtained with 
solutions to the original (a) and rescaled (b) equations, both for the relative residual 10 -7 . 


We next consider a higher frequency, about 3 kHz, corresponding to the wavelength A = 
10 cm. In Fig. 23 we plot convergence histories for original and rescaled equations in this case. 
It is seen that convergence of the rescaled equations is fairly independent of variations - within 
a factor of about two - in the scaling parameter £. 



nn 

Figure 23: Convergence of the iterative solution to the transmission equations with the outer 
head surface, for A = 10 cm and the tissue density p m / p 0 = 1000: no rescaling (£ = 1), and 
rescaling with factors £ = 1/500, £ = 1/1000, and £ = 1/2000. 
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The energy flux distribution computed with rescaled equations is visualized in Fig. 24. 
In this case the solutions to the original equations are close to those shown in the Figures 
only provided the relative residual level ^ 10 -5 or better, while for the rescaled equations the 
relative residual level ^ 10 -3 is entirely adequate. 
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(a) (b) 

Figure 24: Distribution of the energy flux density in a model of a human head filled with a 
homogeneous tissue with p m /p 0 = 1000 and n m — 0.2, and for A = 10 cm, seen from two points 
of view, (a) and (b). 


For comparison, we plot in Fig. 25 the energy flux distribution for an even higher frequency 
(A = 5 cm). It can be seen that it is more concentrated near the entrance to external ear canal, 
even though the ear canal itself is not modeled in the considered geometry; hence, the effect 
should be attributed to the shape of the auricle. 
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Figure 25: The same as Fig. 24, but for the frequency twice as high (A = 5 cm). 


B.5.2 Energy flow in model of a human skull (B) 

Here we show results of computations for a model of an empty human skull. We do not consider 
this problem a realistic simulation of sound conduction in a human head , in which the skull is 
both surrounded by soft tissues and filled with them. This fact has been recognized long ago 
in the experimental work on bone conduction (e.g., [17, 18]); Nevertheless, it is of interest to 
analyze the energy flow in an isolated skull, to compare it with the behavior of a more complete 
model. 

Our skull model is described by a triangular mesh with N v = 112, 038 vertices (hence the 
total number of unknowns is twice as large). The triangulation, shown in Fig. 26, is relatively 
uniform: the minimum, average, and maximum values of the edge length are 0.07 mm, 1.52 mm, 
and 4.62 mm. 

Fig. 27 shows the distribution of the pressure field on the surface of the model, for wave¬ 
length A = lm, corresponding roughly to the frequency 300 Hz. 

As expected for a high-density material (nearly a hard surface), the pressure is of order 
unity. As also expected for a nearly hard surface, the velocity is much smaller (Fig. 28). 

The velocity distribution is seen to be concentrated in the area of the temporal bone, simply 
because the bone is thinner in this region than elsewhere throughout the skull. The energy 
flux, shown in Fig. 29, is also small and concentrated in the same area as the velocity field; 
this is expected, since the distribution of pressure is quite uniform. The distributions of the 
flux density in Figs. 29(a) and Figs. 29(b) show that energy enters the bone from the outside 
and leaves through the inside surface. 
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( a ) 


(b) 


Figure 26: Triangulation of the left half of the skull model, as seen from the interior side (a), 
and a detail of the triangulation in the vicinity of the inner ear (b). 
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Figure 27: Distribution of the absolute value of the pressure field the exterior (a) and interior 
(b) surfaces of the left half of the skull. 
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Figure 28: Distribution of the absolute value of the velocity field the exterior (a) and interior 
(b) surfaces of the skull. 
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Figure 29: Distribution of the energy flux density on the exterior (a) and interior (b) surfaces 
of the skull. Positive/negative values indicate energy entering/leaving the bone. 


For comparison with solution to the high-contrast transmission problem, we show in Fig. 30 
the solution for the pressure field for the hard-surface (Neumann) problem. While the pressure 
on the outer surface (Fig. 30(a)) of the hard-surface skull is comparable to that in the trans¬ 
mission model, it is much (about 10 times) larger in the interior of the skull cavity (Fig. 30(b)). 
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This phenomenon can be intuitively understood by invoking the equivalence of the acoustic 
Neumann problem to the problem of steady-state heat transport in the presence of insulating 
boundaries, in which the pressure becomes the temperature, the Neumann boundary condition 
dp/dn means zero energy flow through the boundary, and the incident plane wave plays the 
role of an external energy source. The “physical intuition” suggests then that the heat should 
be trapped in cavities. 

The above comparison emphasizes the difference between the hard-surface problem and 
a transmission problem, even with a very large density contrast, hence small energy fluxes 
through the interfaces: even a small energy flow through the surfaces prevents its accumulation 
in a cavity. 

As another comparison, Fig. 31 shows the distributions of the velocity field in the soft- 
surface (Dirichlet) problem. It is, of course, very different from the distribution in the trans¬ 
mission problem in Fig. 28, and (in agreement with the expected solution of the equivalent 
electrostatic problem) it exhibits maxima at pointed elements of the surface. 



Figure 30: Distribution of the absolute value of the pressure field on the exterior (a) and 
interior (b) surfaces of the skull modeled as a hard surface (Neumann boundary conditions). 
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Figure 31: Distribution of the absolute value of the velocity field on the exterior (a) and interior 
(b) surfaces of the skull modeled as a soft surface (Dirichlet boundary conditions). 


Finally, Fig. 32 displays some details of the energy flux distribution on surfaces in the area 
of the temporal bone. It can be seen here that the increased energy flow through the outer 
shell of the skull is related to the small bone thickness. One can also see the cochleas (as 
cavities in the bone). As indicated by the Figures 32, and in more detail in Fig. 33, the energy 
flow through the cochlea is directed from the inside to the outside of the skull. 
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Figure 32: Distribution of the energy flux density in the area of the left temporal bone, seen 
from the back and outside (a) and inside (b) of the skull. 
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Figure 33: Distribution of the energy flux density in the area of the left temporal bone, 
including the outer auditory canal, and the middle and inner ear structures. 
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B.5.3 Energy flow in a model of a human head with skull and soft tissues (C) 

We discuss here some results obtained with a model of a human head consisting of a skull 
embedded in undifferentiated soft tissues filling the outer “skin” boundary. In this model the 
head surface includes the outer ear canals terminated with a flat disc-shaped surfaces modeling, 
approximately, the tympanic membranes. A part of the model is visualized in Fig. 10(a). The 
triangulation used in the geometry results in about 260, 000 unknowns. Most of the reported 
computations were done at the frequency of about 6 kHz, corresponding to the wavelength 
A = 5 cm. 

For comparison, we also carried out a number of computations with a head model without 
the outer ear canals. This model can be considered a realization of a perfect ear-plug , with which 
the outer auditory canal is effectively closed. We start with discussing results of computations 
with this model. 

Head surface without the outer ear canal. The outer surface of the head is in this case 
described by the same model as in Figs. 21 and 23. Fig. 34 shows the resulting density of the 
energy flux through the surfaces of the model. The energy flux distribution on the exterior 
(skin) surface is nearly the same as in the model without the skull (Fig. 25(a)). A coronal 
cross-section through the model is shown in Fig. 25(b). The cross-section plane intersects the 
outer auditory meatus and the cochlea (a cavity in the skull bone). Further details of the 
vicinity of the left meatus and the cochlea are shown in Fig. 35, in which the range of the flux 
density has been reduced to [— 0.001, 0.001] (we recall these values are relative to the incident 
flux density). The boundary of the cochlea is visible near the right edge of Fig. 35(a), which 
presents the view from the back and the left of the head. Fig. 35(b) - the view from the right 
- shows positive values of the energy flux, approximately 0.0002, on the right boundary of the 
cochlea. Those values indicate that the energy flows flows from the cochlea (a cavity) into 
the bone. Similarly, the negative flux values on the opposite side of the cochlea show energy 
flowing out of the bone and into the cochlea. The overall energy flow through the cochlea is 
this from left to right. 
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(a) (b) 

Figure 34: Distribution of the energy flux density in a model of head involving soft tissues and 
the skull: energy flux through the exterior surface (a) and through the interior interfaces (b). 
The latter shows a cross-section of the head, seen from the back and right side. 
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Figure 35: Details of the energy flux density distribution of Fig. 34. 


Head surface with the outer ear canal. The following computation were carried out with 
the full head model, including the outer ear canal and, as before, the skull (Fig. 10(a)). 

The significance of the outer auditory canal in air-conducted sound transmission is obvious. 
However, the results given below also indicate that the outer ear canal plays an important role 
in transporting the sound energy into the interior of the skull, and that a relatively large 
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amount of energy transmitted through the tissues emanates from the wall of the canal. Those 
phenomena are illustrated by the pressure, velocity, and energy flux distributions shown in 
Figures 36 and 37. The values of the energy fluxes are now significantly (by a factor of 10 
or more) larger than in the previous model. An interesting feature is that in different parts 
of the canal the energy flows from the air into the tissue or in the opposite direction. Such 
pressure, velocity, and energy flux distributions are characteristic of a standing acoustic wave 
being formed in the ear canal [3]; consequently, these distributions change with the frequency 
(as confirmed by the computations). 

The following Figures 36 and 37 show the distributions of the pressure and velocity fields, 
and the energy flux density in the region of the left temporal bone. 




(a) (b) 

Figure 36: Distributions of the absolute value of the pressure (a) and velocity (b) fields, \p(r)\ 
and |g(r)|, on the material interfaces in the vicinity of the left temporal bone, for A = 5 cm. 
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Figure 37: Distribution of the relative energy flux density, F(r) of Eq. (B.6), on the interfaces 
in the vicinity of the left temporal bone, in a wider (a) and a smaller (b) regions, for A = 5 cm. 


Some features observed in the solution are as follows: 

1. The outer auditory canal supports a standing wave with a node in pressure at about 
half-length of the canal (Fig. 36(a)). By examining the real and imaginary parts of the 
solution one can see that the pressure changes sign at that point. 

2. The pressure p and the velocity field q (i.e., the velocity component normal to the walls 
of the canal, not the longitudinal velocity) has a more complex behavior and is not simply 
related by a proportionality constant (interface impedance). 

3. As a consequence, the direction of the energy flow through the walls of the outer auditory 
canal changes sign along its length. 

As the wavelength increases, the behavior of the pressure and velocity inside the outer 
auditory canal has no longer a wave character. Nevertheless, the fields and the energy flux 
are still largely concentrated inside the canal. The behavior of the resulting energy flux for 
A = 60 cm is illustrated in Fig. 38. 
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Figure 38: Distribution of the relative energy flux density F(r) on the interfaces in the vicinity 
of the left temporal bone for A = 60 cm. 


Energy flow through the cochlea. Based on the solutions for the two head models, we 
compare now the amounts of energy flowing through the inner ear (more specifically, the surface 
of the cavity housing the cochlea). As a quantitative measure of that energy, we consider the 
integral of the absolute value of the average energy flux density over the surface S of the 
cochlear cavity, 



(B.70) 


where 1^1 is the area of the surface. (We recall that, since the energy is conserved the total 
energy flux through a closed surface is zero - a property satisfied also, to a good accuracy, by 
the numerical solution.) 

In computing the flux density (5.2) we use the surface S shown in Fig. 39 or a similar one. 
The surface is open: one of its boundaries opens into the vestibule, and the other to canal 
containing the cochlear nerve (which then passes through the inner meatus); the solution for 
the pressure and velocity fields is, of course, always computed with closed surfaces separating 
various materials. 
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(a) 


(b) 


Figure 39: Discretized surface S of the skull cavity containing the left cochlea, used in the 
computing the energy flux ^ 5 , seen from the front (a) and from the back (b) of the head. The 
average edge length in the discretization is about 1 mm. The dark interior of the surface is the 
outer side of the bone surface. 


Fig. 40 shows the average relative energy flux (5.2) as a function of the wavelength of the 
incident wave. At larger wavelengths the presence of the outer auditory canal increases the 
energy flux by about a factor of two. At smaller wavelengths the effect is larger and exhibits 
a resonant behavior. The first, sharper, peak at A ^ 5 cm is likely to be due to a resonance in 
the outer ear canal. The origin of the broader maximum at A ^ 15 cm is less clear. One might 
speculate that it is related to vibrations of the entire skull or head; but, since it is absent in 
the model without the auditory canal, those vibrations would have to be excited specifically 
by the energy transported to the interior of the head. 
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Figure 40: The average energy flux density through the surface of the cochlear cavity (relative 
to the incident wave flux density), as a function of the wavelength, for the model with the skull 
and soft tissues. 

In order to assess the sensitivity of the energy flux <P S on the details o the inner ear, we 
carried out computations for three different geometries. The results, for the model of the head 
with the ear canals, shown in Fig. 41, indicate that the differences are minor. 



Figure 41: The average energy flux density through the surface of the cochlear cavity (relative 
to the incident wave flux density), as a function of the wavelength, for the model with the skull 
and soft tissues and three geometries (gl, g2, g3) of the inner ear cavity. 


In order to correlate the observed resonant-type of the energy flux with the spatial flux 
distribution, we plot in Fig. 42 the flux density on the surface of the cochlear cavity, for a set of 
selected frequencies. The plots represents the left cochlea, subject to a sound wave incident on 
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the left ear, and seen from the front of the head. The positive and negative flux values indicate 
the energy flowing into and out of the bone surrounding the cochlea. Hence, for example, the 
distribution at A = 3 cm indicates that the energy flows from the right to the left, i.e., from 
the outside to the inside of the head. However, the energy flow changes its direction twice , 
between the wavelengths 3 cm and 4 cm, and again between 12 cm and 14cm. 
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Figure 42: Distributions of the energy flux density on the surface of the left cochlear cavity 
(geometry gl) seen from the front, for a model of a head with the outer ear canals, and for 
the indicated wavelengths A. Ranges of variation of the relative flux density F(r) are shown 
in the square brackets. 
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Figure 43: The same as Fig. 42, but for a more detailed model of the inner ear cavity (ge¬ 
ometry g2), without the skull and with the outer ear canals, filled with the tissue of density 
pjp 0 = 1000. Discretization of the cochlear cavity is somewhat different than in the previous 
computation. 
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Figure 44: The same as Fig. 42, but for an even more detailed model of the inner ear cavity 
(geometry g3), without the skull and with the outer ear canals, filled with the tissue of density 
pjp 0 = 1000. Discretization of the cochlear cavity is somewhat different than in the previous 
computation. 


By looking at the flux distributions on the inner ear cavity and the outer ear canal we 
can verify that the sudden changes in the direction of the energy flow are are due to changes 
of the field distributions in the outer ear canal near resonances seen in Figs. 40 and 41. As 
an example, we plot in Figs. 45 flux distributions on the middle- and inner-ear cavity for 
wavelengths 12 cm (a) and 14 cm (b). 

The above observations are supported by the fact that in the model without the outer 
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auditory canals, and thus in the absence of noticeable resonances, the energy flow is, at all the 
considered wavelengths, directed from the outside to the inside of the skull. 
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Figure 45: Energy flux density distributions on the left middle- and inner-ear cavity for wave¬ 
lengths 12 cm (a) and 14 cm (b), seen from the top. The canal directed to the left and up is 
the inner meatus. 
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Figure 46: Energy flux density distributions on the left outer ear canal for wavelengths 12 cm 
(a) and 14 cm (b), seen from the front. 


B.5.4 The role of acoustic bone conduction: energy flow in the absence of a skull 

(D) 

In order to assess the significance of bone conduction vs. soft-tissue conduction, we also com¬ 
puted the amount of energy flowing through the cochlea in a model of the human head filled 
entirely with a homogeneous material, without the skull structure. In this comparison we 
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considered two models of the outer head surface: without outer ear canals (i.e., with “perfect 
ear-plugs”) and with ear canals. 

Head surface without outer ear canals. We consider here two cases: the head filled with 
the soft tissue (p/po = 1000, n = 0.2) and the bone (p/p 0 = 2000, n = 0.4). We then solve 
the resulting transmission problems with the outer head surface (the “skin”) and a part of 
the middle- and inner-ear ear cavity, including the cochlear cavity. Finally, we compute the 
average energy fluxes (5.2) through the same part S of the cochlea cavity surface as used before 
(Fig. 39). The results are shown in Fig. 47, together with the previous plot of Fig. 40 for the 
model without the outer ear canal, but with the skull. 



Figure 47: The average relative energy flux density through the the cochlear cavity as a 
function of the wavelength, for the model without the skull and two tissue densities ((b) and 
(c)), compared to the result for the reference model (a) with the skull. In this model the head 
surface model has no outer ear canals. 


The relations between the fluxes for larger wavelengths can be plausibly explained: 

- Comparison of the models (a) and (b) indicates that the presence of the higher-density 
skull in the soft tissue may reduce the energy flux due to its transmission through addi¬ 
tional interfaces. 

- Reduction of the flux in the model (c) compared to (b) (by about a factor of two) may 
be simply attributed to the higher density contrast in (c). 

On the other hand, we find it difficult to explain in an intuitive way the behavior of the 
energy flux for wavelengths A < 10 cm, especially for the higher-density model. As an example, 
we computed distributions of the pressure and velocity fields, as well as the energy flux density 
on the surface of the head model (c), for two wavelengths, A = 6 cm and A = 7 cm, for which 
the average flux densities <3? s (Fig. 47) are about 6.5 • 10 -5 and 3.5 • 10 -4 . The pressures for 
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the two wavelengths are quite similar, but the velocities for A = 7 cm are significantly larger 
on most of the head surface than for A = 6 cm. As a result, the flux density for A = 6 cm is 
more concentrated near the center of the outer ear, but for A = 7 cm it is more spread over the 
head surface and overall larger. This fact appears to indicate that the larger value of the flux 
<P S for A = 7 cm is due to energy penetrating through a large area of the head surface, rather 
than through the vicinity of the ear. 

At the same time, the wavelength-dependence of the flux for the model (b) - with the 
density p/p 0 = 1000 - is much more smooth, due to a more smooth dependence of the velocity 
field, but we cannot offer an explanation of that fact. Generally, we find that, as expected for 
the large-contrast problems, the behavior of the pressure distribution is nearly independent 
of the density and relatively weakly dependent on the wavelength; however, the velocity-field 
distribution, and thus the energy flux, vary with both the wavelength and the density in ways 
difficult to predict. 
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Figure 48: Distributions of the pressure and velocity fields, and of the energy flux density on 
the outer surface of head for A = 6 cm (left column) and A = 7 cm (right column), for the 
model without the skull and with p/p 0 — 2000, n = 0.4. 
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While the above discussion pertains to the energy flux on the outer surface of the head, 
the energy flux distributions on the cochlear cavity for several wavelengths and for the two 
considered materials are shown in Figs. 49 and 50. The distributions for the two materials 
are similar and show that, except for the first entries, the energy flux grows monotonically 
with the wavelength and becomes more uniform. We recall that the Figures represents the 
left cochlea, subject to a sound wave incident on the left ear, and seen from the front of the 
head. The positive and negative flux values indicate the energy flowing into and out of the 
bone surrounding the cochlea cavity; hence, the energy flux distributions show that the energy 
is flowing from the inside to the outside of the head - a phenomenon observed also, for some 
wavelengths, in the full model with the skull (Fig. 42). Again, this behavior appears to be 
difficult to explain intuitively without performing the actual computation. 
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Figure 49: Distributions of the energy flux density on the surface of the left cochlea cavity, seen 
from the front, for a model of a head without the skull and with the tissue density p/— 1000. 
The distributions are shown in varying scales for a number of wavelengths, from A = 5 cm to 
A = 80 cm. 
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0.000275 

0.000208 

0.00014 

7.31e-05 

5.85e-06 

-6.14e-05 

-0.000129 

-0.000196 

-0.000263 

-0.000331 

-0.000398 



(e) A = 60 cm 


i 


0.000727 

0.000622 

0.000516 

0.00041 

0.000304 

0.000198 

9.24e-05 

-1.34e-05 

-0.000119 

-0.000225 

-0.000331 

-0.000437 

-0.000543 

-0.000648 



(f) A = 80 cm 


| 


0.000845 

0.000722 

0.000598 

0.000474 

0.00035 

0.000227 

0.000103 

-2.08e-05 

-0.000145 

-0.000268 

-0.000392 

0.000516 

-0.00064 

0.000763 


Figure 50: The same as Fig. 49, but for a the tissue density p/p 0 = 2000. 
The distributions are shown in varying scales. 
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Head surface with outer ear canals. In this problem we also consider the reference full 
head model with the soft tissue and the bone (a), and compare it with skull-less heads filled 
with either of the two materials (models (b) and (c)). The results for the average energy flux 
density on the cochlear cavity are shown in Fig. 51. 



Figure 51: The average relative energy flux density through the the cochlear cavity as a 
function of the wavelength, for the model without the skull and two tissue densities ((b) and 
(c)), compared to the result for the reference model (a) with the skull. Here the head surface 
model contains outer ear canals. 


Again, for A > 10 cm the flux densities in the the three models differ by less than about 
30 % and the relations between them can be explained as in the previous problem. For smaller 
wavelengths we observe a close agreement of the models (a) and (b), but the model (c) continues 
to exhibit more resonance-type structure. 

Now, since the skull-less model (at least for the density p/ p 0 = 1000) yields a flux density 
behavior closely similar to the model with the skull, it is also of interest to analyze the flux 
density distribution on the cochlear cavity - as shown, for the full model, in Fig. 42. The 
results for the present model (with a somewhat different cavity discretization) are visualized 
in Fig. 52. 
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A = 3 cm 

[-0.00011,0.00013] 



A = 4 cm 
[-0.0011,0.0025] 



A = 6 cm 
[-0.0011,0.00062] 



A = 8 cm 

[-0.00083,0.00048] 



A 10 cm 
[-0.0017,0.00095] 



A = 12 cm 
[-0.0060,0.0035] 



A = 14 cm 
[-0.0029, 0.0076] 



[-0.0017,0.0035] 



A = 30 cm 
[-0.0012,0.0024] 



A = 40 cm 
[-0.0014,0.0026] 



A = 60 cm 
[-0.0023,0.0040] 



A = 80 cm 
[-0.0028,0.0046] 


Figure 52: The counterpart of Fig. 42, now for a head model without the skull and with the 
outer ear canals, filled with the tissue of density p/p 0 — 1000. Discretization of the cochlear 
cavity is somewhat different than in the previous computation. 


We observe a remarkable correlation between the distributions in Figs. 42 and 52. In 
both cases the energy flux through the cochlea changes direction in the resonance region. 
Interestingly, however, the directions of the energy flow in the two models are almost exactly 
opposite throughout the entire wavelength range! On the basis of this comparison we conjecture 
that 

- The resonant behavior of the flux at smaller wavelengths is independent of the presence 
of the skull, and is thus due to the outer ear canals. 
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- The overall magnitude of the energy flux inside the model is primarily controlled by the 
outer model surface and is fairly independent of the skull. 

- The skull, however, affects the local distribution of the flux, possibly through reflections 
from the interfaces between the bone and the soft tissue. 

B.5.5 Summary of the results on energy flow 

The main findings of our analysis on energy flow through the human head can be summarized 
as follows: 

1. Transmission problems involving multi-region systems of the types listed above can be 
accurately solved by means of either first- or second-kind equations. Expectedly, second- 
kind equations have a certain advantage in terms of better conditioning and thus a smaller 
number of iterations in the solution. 

2. Solution for the pressure on the exterior air-tissue interface is insensitive to the interior 
regions of the model. Actually, the pressure on that surface can be well approximated 
by the solution of a simpler hard-surface scattering (Neumann) problem. 

3. The value of the normal derivative of the pressure on the exterior surface depends pri¬ 
marily on the density of the tissue adjacent of the surface, but it is also more sensitive 
to the materials and the geometry of the interior of the model. 

4. We have observed some interesting features of the distributions of the pressure and its 
normal-derivative on the exterior high-contrast interface, and the corresponding fea¬ 
tures of the energy flux through the surface. As we discuss below, our problem of a 
high-contrast interface involves features of both the hard-surface (Neumann boundary 
condition) and soft-surface (Dirichlet boundary condition) problems. 

(a) The pressure p tends to be concentrated in depressions and indentations of the 
surface, in particular in canals penetrating inside the object (including the outer 
auditory meatus). Intuitively, this behavior is expected by analogy with the tem¬ 
perature distribution in heat transport in the presence of an insulating boundary - 
formally, a problem identical to acoustic scattering off a hard (Neumann condition) 
surface. 

(b) At the same time, the normal derivative of the pressure, dp/dn , tends to concen¬ 
trate at protruding elements of the surface. Again, this behavior can be intuitively 
understood by analogy with the distribution of the charge density distribution on 
the surface of a conductor, described by the Laplace equation with the Dirichlet 
boundary conditions. 

5. We found that the distributions of p and dp/dn depend in very different ways on the 
surface geometry, hence the behavior of the their product (i.e., the energy flux) may be 
highly nontrivial and difficult to predict without an actual computation. In particular, 
the behavior of the energy flux in the vicinity of the outer ear canal suggests that a 
significant part of the non-airborne sound transmission may come from this area. 
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6. The comparison of the results of models with and without outer ear canals (with the latter 
model representing effects of “perfect ear-plugs”) shows that at relatively high frequencies 
(the wavelengths < 15 cm) a large part of the non-air-borne energy transfer to the inner 
ear is due to the resonant behavior of the wave in the outer ear canals. This contribution 
to the energy flux can be significantly reduced by blocking the ear canals. However, as 
the frequency decreases, the contribution of the energy penetrating through the surfaces 
of the head and the skull becomes dominant; and this energy-transfer mechanism can 
only be suppressed by protecting the entire surface of the head. 
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C Draft of the paper 

“Formulation and Applications of the First and the Second 
Kind Elastodynamics Integral Equations in Simulation of 
Elastic Wave propagation in Human Head 
(part 1)” 


E. Bleszynski, M. Bleszynski, and T. Jaroszewicz 
Monopole Research, Thousand Oaks, CA 91360 

Abstract 

We present formulation and selective applications of the boundary integral method 
formulation in elastodynamics applicable to a general multi-domain object composed of 
piecewise homogeneous material regions characterized by distinct Lame material parame¬ 
ters. The most exterior region may have exclusive acoustic properties. 

The constructed integral equation solver offers a versatile numerical simulation frame¬ 
work providing high accuracy, ease of treating high-contrast interfaces and material proper¬ 
ties, essential in modeling intricate geometrical structures with complex biological material 
properties present in the area of human middle and inner ear. 

C.l Introduction 

We present elements of the formulation and implementation of the elastodynamics integral 
equation solver applicable to a general problem involving an object composed of a number 
of domains characterized by different Lame material parameters. The resulting formulation 
consists of a coupled set of surface integral equations, for two unknown surface vector fields: 
displacement u(r), and the traction t(r) fields. The discrete representation of the integral 
equations is constructed by expanding the displacement u and traction t fields on material 
interfaces in terms of piecewise linear basis functions supported on sets of triangular facets 
sharing a common vertex. We have constructed two different versions of integral equations: 
(a) in the first-kind form and (b) in the second-kind form, for the displacement and traction 
fields expanded in terms of the linear node-based basis functions. We constructed explicit, 
numerically stable expressions for Galerkin matrix elements, of all pertinent kernels appearing 
in both types of integral equations for solid-solid, fluid-solid, and fluid-fluid material interfaces. 

In order to be able to handle the large scale realistic numerical simulations we interfaced the 
solver with a suitable FFT-based matrix compression algorithm which reduces the numerical 
complexity of the iterative solution to N\ogN. The choice of the compression method is 
motivated by its efficiency in treatment of subwavelength problems present at frequencies 
associated with typical external acoustic excitation. The parallelization scheme for the matrix- 
compressed version of the elastodynamics integral equations has been developed as well. 

C.2 Lame differential equations 

We start with the Lame equation for the displacement in elastic medium characterized by 
Lame coefficients A and fi The displacement u(r) and the stress tensor f(r) 

f (r , uS) = A/V r • u(r) + /i[V r 0 u(r, ui) + u(r , cj) 0 V r ] 
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satisfy the following equations 


{(A + /i) V x ® Va; + llV 2 x + pco 2 }u( X, u) = p(x)F(x, lj), 

Vjc • t(x,lo) + p(x)(jU 2 u(x,Uj) = — p(x)F(x, bo), 

The Green function of the Lame equation is the second rank symmetric tensor 

G(R) = C(R)I + V R ® V r D(R) 

C(R) = -gs(R), 

(1 

D(R) = 4o[g s (R) - gc(R)}, 


(ik^ 


with 


R = x — y. 

gi k c R 

4 tvR 5 
Jks R 


gc{R) = 
gs(R ) = 


4t tR’ 

V 2 R g c (R) = -k 2 c g c (R)-5(R), 
V 2 R g s (R) = —k 2 s gs{R) - S(R), 

V 2 r D(R) = -L[k 2 c9c (R) - k 2 g s (R)}, 
-\k 2 c g c {R) = %C{R) - k 2 c D(R ), 

pid r£ g 


V r -G(R) = ^V r9 c(R) = 




A T 2 fl 


^r9c(R), 


We define the stress tensor 


r(y) = A IVy • u(y) + y[V y <g> u(y ) + u{y) 

hj{y) Cijj^ rn dj^u rn ( K y ^), 

Cijkl = + Silbjk), 


y .b 


The two wave-numbers, 


(jj 

k c = -, 


k -- 


(C.l) 


(C.2) 


(C.3) 


(C.4) 


(C.5) 


are related to the longitudinal (compressional) and transverse (shear) wave speeds as follows 

(C.6) 


‘\ + 2p 


C S “ 
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f(x,y) = n(y) ■ E(R), 
d 

r%j{x,y) — — Gip{x y)Cjk pq n k (y) = \ri(y ) • S(x y)}ij = ^i{y > )Gij rr iid m Gi k ( < R')^ 

OVq 

ZJ(R)} = A IV • G{R) + /iV ® G(R) + /xG(R) ® V, 

[A7(R)]^-fc = A 5ijd rn G mk (R) + + djGik(R)] = Cij rn id m Gi k (R), 

r jk(x,y) = ni{y)S ijk , 

$ jk (x,y) Sij k Tii{x s ) 

£(R) = -£(-R). 


The Green tensor functions G(R) and £(R) satisfy the following equations 

[(A + /i)V H ®V R + + pu 2 ]G{R) = - JCR)/, 

• £(J2) + pu 2 G{R) = -<S(J2)J, 

T(R) • + ^ 2 G(R) - A(V fi • -Vfl®)V fl • G(R) = -S(R)I , 

d4Qj&AG/ m (/£)] + pcu 2 Gj m (R) = — 8j m 5{R ), 

+ pcj 2 Gj m (R) = — Sj rn S(R ), 

V* •<?(«) = XT2^ Vh5c(jR) ’ 

while t(x) satisfies the following equation 

p(r)u 2 u(r) = 0 , 

^[V r (E> w(r) + «(r) <g> V r ] - f w (r) = 0, 


(C.8) 


(C.9) 


where 


Cijkl = ^ij^kl + H{bik&jl + 6i/5jfc) 


Dijkl ~ 4//, 

1 

4/i 




2A 


- y &, - 3A+2fi 


3A + 2^‘ iSu 

SijS k i 


(c.io) 


Cijkl Cjiki Ciji k C k Hj. 

is the 4-th rank elastic stiffness tensor. A useful relation is 


C , ijifez5jfeTxz(®)5 i G J - m (R) - C ijk id i u j (x)d k Gi rn (R ) = 0. (C.ll) 


C.3 Representation formulae in elastodynamics 
C.3.1 Representation formulae for the displacement field 

We consider a single region i? bounded by dfl containing homogeneous material with Lame 
parameters A and p (see Fig. 53. By assuming that, inside L? the displacement u(x) satisfies 
the Lame equation (C.14) 
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Figure 53: A region f? = fi m separated by df2 in i?o used in the derivation of the integral 
representation for u. 


The integral representations for the displacement and traction fields inside and a 
region can be derived by applying the of Gauss divergence formula to the equations satisfied 
by the Green tensor functions G(R) and Z’(ii) for a single frequency uo\ 

Vy • f(y) + pu 2 u(y) = —p(y)F(y), 

V y • E(y, x ) + pu 2 G(y , x) = -8{y - x )/, 

By multiplying the first of the above equations by G(x,y) from the right side and by 
multiplying the second by u(x) from the left side and by subtracting them we obtain 

Vy • t (y)] • G(y, x) - u(y) • [Vy • E(y, x)} = u(y)S(y - x) - p(y)[F(y) ■ G(y, x)}, (C.13) 

But 

Vy • U(y) • G(y, *)} - u(y) ■ E(y, x)]} = [V y • t( y)\ ■ G(y, x) - u(y) • [Vy • E(y, x)]+ 
+[r(y)] • -[Vy ® G(y, x)\ - [Vy <g> u{y)\ ■ \E(y, *)] 

(C.14) 

where' denotes contraction in the first two indices, 

Also, 

E(x -y) = ~E(y - x), 

G(x,y)\ = G(y,x), (C.15) 

Vy • G(y,x)] = V x ■ G(y,x), 

and 

[Vy ® u(y)} • \E{y,x,u)] - [r{y)\ ■ \V y ® G(y,x )] = 0 (C.16) 

Let us verify the above identity in the indicial notation 

duj ^ ^ . , ^dG jk (y,x) _ duj^ dG nk (y,x ) ^ du n dG jk (y,y) _ 

dyi ijk[V,X) Tij[V) d Vi d Vi ijmn dy m ijmn dy m % 

_du nri dGj k (y, x) du n dG jk (y, x) _ n 

— 0 mnij 0 ijmn ^ ^ — U 

dy m d yi dy m % 

(C.17) 
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We may now rewrite (C.14) as follows 

= • (f(y) • G(x,y)+u(y) ■ £(x,y)]} = u(y)5(x - y) - p{y)[F{y) ■ G(x,y)\, (C.18) 


By taking the volume integral of the above equation we obtain and by applying the Gauss 
divergence formula we may convert left side to the surface integral 

J dW y ■ {t (y) • G(x, y) + u(y ) • £(x, y)j = 
n 

= / dS v {[h(y) ■ f (y)] • G(x, y) + [u(y) ■ n(y) ■ S(x, y)]} = (c. 19 ) 

on 

= [ dS y {t(y)-G(x,y) + u(y)-P(x,y)} 
dQ 

We obtain the integral representation for the displacement field for a region h? bounded by dfl 
which allows to to find the value of the displacement field at any point x in Q in terms of the 
integral of 2 fields u(x) and t(x) on the region boundary for u(x) 


[ dS y [t(y) ■ G(x, y) + u(y ) • P(x, y)] 

an 

t(y ) =n(y) ■ f(y), 

P(x,y) = n(y) ■ £(x,y), 


u(x) for x ^ 17, 
0 for x G 


(C.20) 


The symbols appearing in the integrand of the above integral representations are as follows 

• u(y) is the displacement vector field 

• t(y) is the traction vector field related to the stress tensor r(y) as follows 

i(y) = n(y) • t (y) = 

= An(y)[V y • w(y)] + y{[n(y) • V y ]w(y) + V y [n(y) • w(y)]}, 


• G(cc — y) is a second rank symmetric tensor, the Green function of the Lame equation 
given by (C.2) 

• r(x,y) is a second rank non-symmetric tensor, related to 

• E(x), a third rank stress tensor (symmetric in first 2 indices) 


C.3.2 Representation formulae for the traction field 

In this section we construct the representation formulas for the traction field which we will 
use in the context of first kind of integral integral equation formulation. We consider a single 
region Q bounded by dQ By assuming that, inside Q the displacement u(x) satisfies the Lame 
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equation (C.14) we may construct the integral representations for the traction fields 
inside and a region i?. We obtain the integral representation for the displacement field for a 
region Q bounded by dQ which allows to to find the value of the displacement field at any 
point x in Q in terms of the integral of 2 fields u(x) and t(x) on the region boundary 

f d Sv [Kv) ■ #(*> y) + u (y ) • w (*> y)} = 

dQ 

t{y) = h(y) ■ r(y), 
f(x,y) = h{y) ■ E(x,y), 

where i?+ denotes the position locate d just above the surface Q. 

In order to construct the above integral representation for the traction field outside 
region Q bounded by dfi we applying the traction field operator to the integral representation 
for displacement u(x) 

t[u{x)] = A/[V® • u{x)\ + fi[V x 0 u(x) + u(x) 0 V*] = 

= t{ f dSy [u(y) ■ f ( x, y) + t(y) ■ G(x, y)\ } 
dQ 

= [ dS y [r[u(y) ■ f(x,y)\ +r[t(y) ■ G(x,y)]\ 

dQ 

f[u(y) • f(x, y)} = XIV x ■ [u(y) ■ f(x, y)} + yV x <8> [u(y) ■ r(x, y)] + y[u(y) ■ r(x, y)] <g> V* 

r[t(y) ■ G(x, y)\ = XIV x ■ [it(y) ■ G(x, y)\ + yV x ® [t(y) ■ G(x, y)\ + y[t(y) ■ G(x, y)] ® V* 

(C.23) 

The contributions to traction field components associated with the last 2 terms of (C.57) are 
h(x) ■ f[u{y) ■ f(x, y)] = u(y) • W (x, y) = 

= An(*){V £C • [u(y) ■ f(x, y)]} + y[h(x) ■ V x ][u(y) ■ f(x,y)] + yV x {u{x) ■ f(x,y) ■ h{y)} 
h(x) ■ f[t(y) ■ G(x, y)] = t(y) ■ S(x, y) = 

= Xh(x){V x ■ [t(y) ■ G(x,y)]} + y[h(x) ■ V x }[t(y) ■ G(x,y)\ + yV x {n(x) ■ G(x,y) ■ t(y)} = 
= t{y) ■ {AfVa, • G(x, y)\ <g> h(x) + y[n(x) • V x ]G(x, y) + y[h{x) ■ G(x, y)\ ® V*} 

(C.24) 

The above equations can be used to identify the Green function operators W(x,y) and 
S(x,y) 

W(x, y) = A [P(x, y) ■ V x ] ®h a + y[h a ■ V x \P(x, y) + y[P(x, y) • n a )] <g> V* 

S(x, y) = A^* • G(x, y)\ <g> n(x) + y[n{x) • V x ]G(x, y) + y[n(x) ■ G(x, y)\ <g> V x 

and their matrix elements. 


t(x) for x ^ 17, 
0 for x E 


( r , 99^ 
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We may deduce the form of W (a?, y) from as follows: 

W(x,y) = A [f(x,y) ■ V x ] 0 n a +/i(n a • V x )f(x,y ) + y[f(x, y) ■ n a )] 0X7^ = 

A + 2 ® ^’ a SC (R) "h ^ ^ { fi'b ' ^ x)^x ® Ra d" (fia ' ^£c)^6 ® j - 5 , C'(-^)d' 

+ ^ 2 {(n a • V x )(h b ■ V x )G(x,y) + (n a • V a! )V a: 0 [n 6 • G(*,y)]+ 

+ ( 1 U 5 Vj) [ri a • G(x, i/)]®Vj; + [ri;, • G(x, y ) • fi a ] V *8 VjJ 

(C.26) 


C.3.3 The acoustics limit 

As a check of the representation formulae and the Green functions in elastodynamics, we 
consider in the following their acoustics limit, y —> 0 . 


The Green functions in the acoustics limit. In the acoustics limit the Fourier transform 
of G becomes 


Gij (l-i—tO) 


r --(* qiq h 1 1 1 

lJ y [dij kl ] q>-kl + yklql-kl 

1 




QiQj 


A h*™ q 

Qk 


2 _ b.2 > 

C 


^(R)ijk(n^O) $ij o u2 

q ~ K c 


G(R )„=0 = ~jp i is3 ( R ) + Vfi ® V R g c m, 
c 

S{R) = A IV R ■ G(R) = IV R gc(R), 
f(R) = n(y) 0 V R g c (R), 

S(R) = 'V R g c (R) 0 n(x), 

W(x, y) = -Xk 2 c h b 0 h a g c (R). 


(C.27) 


(C.28) 


Representation formulae in the acoustics limit. We consider a single region 17 bounded 
by dl7 By assuming that, inside 17 the displacement u(x) satisfies the Lame equation (C.14) we 
may construct the following integral representations for the displacement and traction 
fields in a region 17 bounded by dl7 which allows to to find the value of the displacement field 
at any point x in 17 in terms of the integral of 2 fields u() and t(x) on the region boundary 

u(x) = 

dn 

t(x)= j dSy [u(y) ■ W(x,y) +t(y) ■ &(x,y)\ 
an 


JdS y [u(y) 


R(x, y) + t(y) ■ G(x,y) 


(C.29) 
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Since the operators appearing in the integral equations are obtained as boundary limits of 
the expressions for the fields, the “contact” 6(R) terms do not contribute to the expressions 
for the operators. After discarding such terms, we obtain 

G(x,y)^ =0 = - —£ 2 - [Vk<8> Vngc(R)], 

c 

r{x, y ) M=0 = h{y) 0 V Rgc(R), 

y ) M =0 = Vr9c(R) 0 n(x), 

W(x, y )^ = o = -A k 2 c h{y) 0 n(*)0c(i2), 

1 (0.31) 

(bGa )^ =0 = ■ 6)(V« • a)sc(i2), 

C 

(bfa)^ = (b ■ fib)(a ■ V R )g c (R), 

(bSa )^ 0 = (a • n a )(b ■ V R )g c (R)], 

(bWa)^ =Q = -A k 2 c {n b ■ b)(n a ■ a)g c (R) 

C.4 Construction of coupled surface integral equations for piecewise homo¬ 
geneous media for displacement and surface traction fields 

Surface integral equations - or boundary integral equations, BIEs - are applicable to piecewise 
homogeneous materials, and provide solutions for the displacement and traction fields defined 
on interfaces separating different material regions. Fields in the individual regions are described 
in terms of the appropriate Green functions for elastic materials. 

C.4.1 Structure of a system of regions and the resultant integral equations 

As an example, we give below explicit expressions for a system of surface integral equations 
describing a set of homogeneous regions f2 m separated by interfaces; 



Figure 54: A schematic representation of regions i? and interfaces S appearing in surface 
integral equations . 
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Figure 55: A simplified schematics of the topological structure of regions and interfaces . 



Figure 56: A simplified schematics of the topological structure of regions and interfaces in head 
model. 


One of these regions, i?o, is the unbounded background medium (air). The displacement 
and traction fields are assumed to be continuous across the interfaces. 

The resulting system of integral equations is then obtained in the following three steps 

(a) using a conventional representation theorem (found e.g. in Morse and Feshbach) for each 
material region i? for the displacement field in this region written in terms of the surface 
integrals over dQ with “surface sources” displacements and traction fields 

u(x) = + / dS y [u(y) ■ f m (x,y) + t(y) ■ G m (x,y),] 

an 

for the displacement u(r) and the traction t{r) = n(r) • f (r) fields 

(b) imposing boundary conditions on the continuity of u(r) and and of t(r) on each inter¬ 
face (oriented surface) S mn separating the regions two equations per interface (oriented 
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surface) S mn separating the regions (on the negative side of the interface) and Q n 
(on its positive side), 

U m (r) = Un (r), 

mV ) h C.32 

tm(r) = t n (r). 

(c) by writing two two suitable equations following from the boundary conditions for each 
interface S m n • 

C.4.2 Boundary conditions 

If the n-th interface is an interface between 

• two elastic media (1) and (2), 

• two viscous fluids, or, 

• a viscous fluid (1) and an elastic medium (2), 

we impose the following boundary conditions (six in general) 

• equality of the three components of the displacement field on both sides of the interface, 

u^\r) — u^ 2 \r ), (C.33) 

• equality of the three components of the traction field on both sides of the interface, 

r (1 )«W(r) = f (2) w (2) (r). (C.34) 

On the interface between non-viscous fluid and the elastic medium we impose the fol¬ 
lowing reduced number of boundary conditions (four in general): 

• equality of the normal component of the displacement field on both sides of the interface, 

n^\r) • u^\r) = n^(r) • u^ 2 \r ), (C.35) 

• equality of the normal component of the traction field on both sides of the interface. 

n^(r) • t^\r) = n^(r) • t^ 2 \r ), (C.36) 

• vanishing tangential components of the traction field 

n( x )(r) x [n^ 1 ^(r)x)t^ 1 ^(r)] = h^ 2 \r) x [n^ 2 ^(r)x)t^ 2 ^(r)] = 0, (C.37) 

on both sides of the interface 

Finally on the interface between two non-viscous fluids we impose the following reduced 
number of boundary conditions (four in general): 
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• equality of the normal component of the displacement field on both sides of the interface, 

n^(r) • u^\r) = n^\r) • u^ 2 \r ), (C.38) 

• equality of the normal component of the traction field on both sides of the interface. 

h^ l \r) • t^\r) = n^(r) • t^ 2 \r ), (C.39) 

• vanishing tangential components of the traction field 

n( x )(r) x [n^ 1 ^(r)x)t^ 1 )(r)] = h^ 2 \r) x [n^(r)x)t^ 2 ^(r)] = 0, (C.40) 

on both sides of the interface 

C.4.3 Choice of the integral equation formulation 

The choice of the integral equations is not unique. It is possible to write several different set of 
integral equations by taking different linear combinations of equations representing boundary 
condition for continuity of displacement and and traction fields across material interfaces. 
While all such integral equations are theoretically equivalent, they tend to differ in terms of 
accuracy, computational resources needed an solution convergence 

We considered the two integral representations for the displacement and of the traction 
fields 


u(x) = f dSy[u(y) ■ f(x,y) + t(y) ■ G(x,y)\ 
da 

t(x)= j dS y [u(y) ■ W(x,y) + t(y) ■ &(x,y)\ 
da 


(C.41) 


C.4.4 Construction of the second-kind integral equations system 

On the interface S mn we write 6 equations for 6 unknowns u(x),t(x) 

lim u(x) in regioni? m 
xedQm 

lim u(x) in region!?^ 

xedftn 

The resulting second kind pair surface integral equations for u and t for the 
interface separating two regions fi m and Q n is 
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\u(x)+ J dSy [u(y) ■ f m (x, y) + t(y) ■ G m {x, y)\ 

S 

mn 

~ f dS y [u(y)-f m (x,y) + t(y)-G m (x,y)]=5 m0 u m (x) 

S . (zdf2 ry 
tm m o . 
i^n irn 

\ u(x) - J dSy \u(y) ■ f n (x, y) + t(y) ■ G n (x, y)\ 
s 

mn 

+ ^2 - J dS y [u(y) ■ f n (x,y)+ t(y) ■ G n (x,y)\ = 5 n0 u m (x) for x € S n 


S . £df2 rt 
nj n o 


(C.42) 


C.4.5 Construction of the first-kind integral equations 

The resulting first kind pair surface integral equations for u and t 

On the interface S mn we write 6 equations for 6 unknowns u(x),t(x) 

u m {x) = u n (x) for x on interface S mn 
t m (x) — t n (x) for x on interface S mn 

The equations of our interest can be obtained in the following way: 

By taking the difference of the equations (C.45) we obtain the first equation of the second 


set 


J dS y {u(y) • [f m (x,y) + f n (x,y )] +t(y) ■ [G m (x,y) + G n (x,y)]} 

S 

mn 

~ ^2 J dS y {u(y) ■ [f m (x, y) + f n (x, y)] + t(y) ■ [G(x, y) + G n (x, y)]} = (c.43) 

1 

[^mO - 5 n 0 ] «“(*) 


S . £ d f2 

im m o . 
i^n irn 


By applying the traction filed operator n( x)) ■ {t{u{x)\ to both sides of the above equation 
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(C.44) we obtain the second equation of the second set 


J dS y {u(y) ■ [W m (x,y) + W n (x,y)] + t(y) ■ [& m (x,y) + & n (x,y)\] 

S 

mn 

- / dS y{ u (y) ' \Wm{x,y) + W n (x,y)} +t(y)- [S(x,y) +S n (x,y)]} (c 44) 

“ si 

= PmO - 5 nol *“(*) 


S ■ (zdf2 ry 

zm m S . 
i^n irn 


The resulting first kind pair surface integral equations for u and t for the interface 
separating two regions and Q n is 


\u(x)+ J dSy \u(y) • r m (x, y) + t(y) ■ G y)] 

S 

mn 

~ J dS y [u(y)-f m (x,y) + t(y)-G m (x,y)]=5 m0 u m (x) 

S. edf2 c 
zm m S . 
z^n im 

\ t(x) - J dSy [u(y) ■ W n (x , y) + t(y) ■ S n (x, y)] 
s 

mn 

+ - J dS y [u(y) ■ W n (x,y) + t(y) -<P n {x,y)] = S nQ t m {x) for x G S mn . 

S ■ G df2 ry 
nj n o ■ 

• / ni 

jytm J 

(C.45) 

With reference to Fig. 55, the first of the above integral equations represents contributions 
to the displacement field u on the interface S mn due to the displacement and traction fields u 
and t on the same interface (the first integral) and on other interfaces, S^ m , forming boundaries 
of the region Q rn with other regions 17^, z ^ n. 

The integrals involve Green functions 


Gi,ri,Wi, and z = l,...,7V 


describing propagation of the displacement and stress fields in the region fij. 

Similarly, the second of the above integral equations represents contributions to the field u 
on the interface S mn due to the fields on the boundaries of the other region, i? n , adjacent to 
the interface. The r.h.s.s of the above equations are the incident fields due to distant sources 
in the region i?o (hence the delta-functions 8 m o and 5 n o)- 

The details of the construction and regularization of T), 4?^, and^ are described in the 
subsequent section. 

We discuss here briefly the form of the surface boundary equations and their discretization, 
which is now being implemented in our solver. 
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C.5 Basis functions and discretization of surface integral equations 

In order to solve the surface integral equations numerically, it is necessary to make assumptions 
on the discretization of the solution, i.e., on the trial basis functions, and on the test basis 
functions. 

In our implementation we use a discretization uniquely determined by our choice of dis¬ 
cretization in the surface integral equations We also use, similarly to the volumetric problem, 
the Galerkin discretization, i.e., identical trial and testing basis functions. 



Figure 57: A set of six triangles sharing a common vertex v a supporting the node based basis 
function 


In the surface problems we assume the displacement field is expanded in piecewise linear 
basis functions supported on sets of triangles. The resulting surface basis functions are piece- 
wise linear vector basis functions describing the components of the displacement field u. By 
symmetry between the displacement and traction fields in the integral equations, we assume 
analogous linear basis functions for the components of t. 

According to the above criteria, we specify the basis functions as follows: 

For each vertex v a of the surface mesh we define three vector basis functions, denoted 
/ 0 a (r), representing displacements in the x, y , and z directions. Correspondingly, the index a 
refers to the vertex and the direction, a = (v a , m), m = 1, 2, 3 (or m = x, y, z). 

Each such function, ^ a (r), is associated with a vertex v a and supported on a set of triangles 
(facets) f a sharing that vertex. We parameterize the basis function as 

^ cX r ) = ^ v a ,m ( r ) = e m ( t>v a ( r ) > (C- 46 ) 

where e m is the unit vector along the m-the axis, and is a scalar basis function defined by 

0a(*) = O) = J2 <,,/«(*) ’ ( C - 47 ) 

fOi^J~ a. 

where the sum is taken over the set T a of all facets f a sharing the vertex v a . Further, each 
of the linear functions (j) v ^ supported on the facet / a , is uniquely defined by setting 
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its value to unity at r — v a and to zero at the remaining vertices of the facet, 
expression is 


K,f a ( x ) = 


i - 


V a ■> fa 


n 


Va 5 fa 


(x-v c 


*/«(*) 


An explicit 


(C.48) 


where Xf a ( x ) is the characteristic function of the facet / a , h v ^ ^ is the unit outer normal 
to the facet edge opposite the vertex v a , and h Va j a is the facet height relative to that edge. 
Components of the basis function are then 


C(0 = = S mi (f> Va {r ) 


(C.49) 


In order to simplify the notation in the following we will use a simplified notation 


A( x ) = ^v a ,m{ x ) = AiAJx) = M x )Xa(x), 

a = v a ,m 


(C.50) 


We represent the linear basis functions components associated with one triangle a or b as 
follows: 


^ dipt) — 

^b{ x) = b<f> b (x), 


<t>a{x) = ^aXa{x) = 

4>b{x) = ibXb{x) = 


1 - 


1 - 


h a -(x- Va) 
h-a 

hb • {x - Vb ) 
h b 


Xa(x), 

Xb{x). 


(C.51) 


As follows from the construction, the scalar and vectorial basis functions (C.47) and (C.46) 
are two-dimensional analogues (actually, restrictions) of the piecewise linear basis functions 
supported on tetrahedra and used in our previously developed volumetric formulation The 
advantage of this discretization scheme is that the solutions of the surface and volumetric 
equations can be directly compared with one another. 

Let r ) and 'ipp(r ) be the two vector basis functions associated with nodes v a and y#, 


= ot<j> Va (r), 


y3<V) = /3<MO- 


a = 


0 = 


a y 

Oi z 

fix 

0y 

Pz 


(C.52) 
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C.5.1 Matrix elements of the integral equation kernels appearing in representa¬ 
tion formulas for the displacement and traction fields 

The explicit forms of the matrix elements of G(x,y) , r(x,j/), and W(x,y) sand¬ 

wiched between single triangle contributions to the linear basis functions ^ a (x) and 3q,(y) are 
: associated with triangles a and b 

&a{ X) = a<f> a (x), 

&b{x) = b<f> b (x), 

M x ) = CaXa(x) = [1 - — ^ - —}Xa(x), 

M x ) = &x&(*) = [i - — ^ — }xb{x). 

are 

[&b(y) ■ G(x, y) ■ # a (x)\ = <f> b (y)<f> a (x)[b ■ G(x, y) ■ a] 

[&b(y) ■ r(x, y) ■ #a{x)] = cf>b(y)(f>a(x)[b ■ f(x, y) ■ a] 

[&b{y) ■ &(x, y) ■ #a{x)\ = (j) h (y)(j) a (x)[b ■ ${x, y) ■ a] 

\&b(y) -W(x,y) ■ V a {x )] = 4> b {y)4> a {x)[b ■W(x,y ) • a], 

with 

[b - f ■ a} = \(b - h b )(a ■ G ■ V x ) + y{h b • V a .)(o • G ■ b) + y(b ■ V x )(a ■ G ■ n b ) 

[b-S - a}= A (n a ■ a)(b ■ G ■ V x ) + y(n a ■ V R )(a ■ G ■ b) + y(a ■ V n)(n a Gb) 

[b-W ■ a] = X[b ■ f ■ V^Ka • h a ) + y(n a ■ V x )(b ■ f ■ a) + /i(a • V x )(b ■ f ■ n a ) 

\ 2 k 2 

= ~ A + 2/j ( a ' ha ^ b ' ™b)gc{R)+ 

+ • a)(h b ■ V x )(b ■ V.) + • 6)(n a • V I )(a • Va.)]5c(-R)+ 

+ /i 2 {(ri a • V x )(fi b ■ Va:)[6 • G(x,y) ■ a] + (h a • V x )(b ■ V x )[rift • G(x,y) ■ a} + 

+ (a ■ V x )(h b • V x )[b ■ G(x, y) ■ n a ] + (a • Va,)(6 • V x )[h b ■ G(x,y) • n a ]} 

(C.55) 


(C.53) 


(C.54) 
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By substituting explicit expressions for G in terms of gs and gc we obtain 


(bGa) = (a ■ + (a ■ V*)(6 ■ V r ) Ss{R) T/ c(B) , 

fi [ikg 

(bfa) = - — (n 6 • 6)(a • Vr)3c(-R) + {(a • &)(n& • Vr) + (n fc • o)(6 • V R )} 3 5 (i?)+ 

+ Tj( a ' ^.R)(& • V.R)(n& • Vr)[ps(.R) - 5c(-R)], 

K s 

(bSa) = — + (»a • a)(6 • Vr) 3 c(-R) + {(a • &)(n a • V fi ) + [n a ■ 6)(a • Vr)} 5 s(.R)+ 

+ 7^( a ‘ Vr)(& ‘ ‘ VR)[5s(il) - 5 c(-R)], 

fc 5 

\2h,2 

(b-W-a) = - - - ° (a • n a )(b ■ h b )g c {R)+ 

^ _|_ 2^ [(^ a ' ^)ix^b ‘ V x )(b ■ V x ) + • 6)(n a • Vj)(o • V£ C )]5'c(-R)+ 

+ /x{(a • 6)(n a • Va.)(n fe • V x ) + (n fe • o)(n a • V a ,)(6 • V x )+ 

“I - (^a ' b)(<l ■ Vx){Rb V x ) ~t~ (^a ' R&)(® V x )(b • V x)}fl , S'(7?) + 

+ I(«a • v x )(n 6 • V x )(a • V*)^ • V*)[0s(i2) - <7c(#)]- 

k g 

(C.56) 

C.5.2 Construction of the integral representation and matrix elements for the 
traction field 

We consider the second choices of coupled integral equations for a multi-region problem. An 
apparent advantage of this formulation is that it enforces boundary conditions on each interface 
in a stronger way than the first set. 

In order to construct the integral representation for the traction field in a region 
Q bounded by dQ we applying the traction field operator yo the integral representation for 
displacement u(x) 

t[u{x)\ = A /[Vx • u(x)\ + /i[V* ® u(x) + u(x) <g> V*] = 

= t{ f dSy \u(y) ■ f(x, y ) + t(y) ■ G(x, y)\ } 
dn 

= f dS y [t[u(y) ■ f(x,y)]+t[t(y)-G(x,y)]] 

dQ 

r[u(y) ■ f(x, y)\ = A IV x ■ [u(y) ■ f(x, y)\ + yV x ® [u(y) ■ r(x, y)\ + y[u(y) ■ r(x, y)\ <g> V x 

r[t(y) ■ G(x, y)\ = XIV x ■ [t(y) ■ G{x, y)\ + yV x ® [t(y) ■ G{x, y)\ + y[t(y) ■ G(x, y)] <g> V x 

(C.57) 
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The contributions to traction field components associated with the last 2 terms of (C.57) are 
h(x) ■ f[u(y ) • f (x, y)\ = u(y) ■ W ( x , y) = 

= Xn(x){V x ■ [u(y) ■ f(x,y)]} + y[n(x) ■ V x ][u(y) ■ f(x,y)\ + yW x {u(x) ■ f(x,y) ■ n(y)} 
h{x) ■ f[t(y) ■ G{x, y)\ = t(y) ■ S(x, y) = 

= Xh{x){V x • [t(y) ■ G(x,y)]} + y[h{x) ■ V x ][t(y) ■ G(x,y)\ + yV X {n(x) ■ G(x,y) ■ t(y)} = 

= t(y) ■ {A[V X • G(x, y)\ 0 n( x) + y[h(x) ■ V x ]G(x, y) + y[fi(x) ■ G(x, y)\ 0 V*} 

(C.58) 

The above equations can be used to identify the Green function operators W(x,y) and 
S(x,y) 

W(x, y) = A [f(x, y) • VJ 0 h a + y[h a ■ V x ]f(x, y) + y[f(x, y) • n a )] 0 V x ^ 

S(x, y) = A[V X • G(x, y)\ 0 h(x) + y[h{x) ■ V x ]G(x, y) + y\fi(x) ■ G{x, y)\ 0 V x 
and their matrix elements. 

The explicit forms of the matrix elements of W(x,y) and ^(x,y) G(x,y) sandwiched 
between ^ a (x) and 'I'biy) are : 

[&b(y) -w(x,y) • Shj(sc)] = <t>b(y)<t>a(x)[b-W(x,y) ■ a} = {n(*) • t[& b (y) ■ f(x,y)]} ■ & a {x) = 
= cf>a(x)(f>b(y){X(a ■ n a ){'V x ■ [b ■ f(x,y)}}+ 

+ n[h a ■ V x ]{[6 • f(x, y)] ■ a} + y(a • V x ){[6 • f(x, y)] ■ h a )} 

l^b(y) -&(x,y) ■ &a(x )] = (f) h (y)(f) a (x)[b ■ S(x,y) ■ a] = {n(x) ■ f[^ b {y) ■ G(x,y)}} ■ & a (x) = 

= cf> a (x)(f> b (y){\(a ■ n a ){b ■ [V x • G(x,y)]} + y[h a ■ Vja • G(x,y) ■ &]+ 

+ y(a • V x ){n a • G(x, y) ■ b)} 

(C.60) 

We may deduce the form of W(x, y) from as follows: 

W(x, y) = A [f(x, y) ■ V x ] 0 h a + y,(h a ■ V x )f(x, y) + /jl[P(x, y) ■ n a )] 0 V x = 

AA: 2 2A u 

— A + 2(i^ lh ^ ^ a 9c{-^) T ^ _|_ ^ji ^' ^ £C )^ £C ® ^ x)^b ® ^ x}9 c(R )T 

+ y 2 {(h a ■ V x )(n b ■ V x )G(x,y) + (n a ■ V X )V X 0 [n b ■ G(x,y)} + 

+ (n b ■ V x )[n a • G(x,y)\ 0V* + [n b ■ G(x,y) ■ h a ]V x 0 V x } 

(C.61) 

By using the following auxiliary expressions for terms appearing in the matrix elements 
' G(x, y) — — — 2 Vx9c{R)i 

v x • [V* • G(x,y)} = [V* • G(x,y) ■ V x ] = -J^-\g c {R) + S(R)}, (C.62) 

V x • [G(x, y) b] = [b- G(x, y) ■ V x ] = ■ V x )g c (R), 
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we find the following formula for the matrix element of W 


b • W(x,y) a = 

a 2 P 

= ~ A + 2 M ( a • ^)( b • nb)gc{R)+ 

+ ^ [(n a ■ a)(n b ■ V x )(b • V x ) + (fi b ■ b){h a • V £C )(a • V x )]gc(R)+ 

+ H 2 {(h a ■ V x )(h b ■ V x )[b ■ G(x,y) ■ a] + (n a ■ V x )(b ■ • G(x,y) ■ a] + 

+ (a • V x )(n b • Va.)[6- G(x,y) ■ n a \ + (a • V x )(b- • G(x,y ) • n a ]} 


b -W(x,y) ■ a = 

\ 2u2 

~ x + 2 ~ ( a ' n a )(b ■ n b )g c (R)+ 

+ -^p^-[(w a • a)(fi b ■ V x )(b • V.) + (n b ■ b)(n a ■ V x )(a ■ V x )]gc(R)+ 

+ M{(« ' b )(n a ■ V x )(n fe • V*) + (n 6 • a)(n a ■ V x )(b ■ V x )+ 

+ (n a • b)(a ■ V x )(h b ■ V x ) + (n a • n b )(a ■ V x )(b ■ V x )}gs(R)+ 

d" T2"(^a ’ ^ x){f^b ' ’ ^x)(b ’ ^x)[9s(R) ~ 9c(R)]- 

k S 


C.5.3 Basic formulae needed for the evaluation of matrix elements of G, T, Z, 
and W 

The kernels G, T, and W need to be sandwiched between the piecewise linear node 
basis function associated with a set of triangles sharing common vertices a and (3. 

We consider a pair of triangles a and b with normal unit vectors fi a and fib The normal 
vectors to the sides a and b are h a and hb We note that n a is perpendicular to the edge l a 
facing vertex a. 

Let us represent the linear basis functions components associated with one triangle a or b 
as follows 


&a{x) = a<f) a (x), 
& b { x) = b(f> b (x), 

(paix) = £aXa(x) 
M x ) = &Xb(x) 


1 - 


hg^jxjz V °) 

K 

hhj_{xj- v b ) 

h b 


Xa(x), 

Xb(x). 


(C.65) 
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We may decompose both a and b into tangential and normal components to triangle faces 

a = (h a x a) xh a + (n 0 ■ a)h a = a t + a n 
b= (h b xb) xh b + (h b ■ b)h b = b t + b n 

Q”ii = ' a)fl a 

at = {n a x a) xh a = a t a t (C.66) 

= ("aa ' a) 
a t = \/l a | 2 - On) 
l a = h a x h a 


Tangential derivatives of basis functions (we neglect the contributions associated with deriva¬ 
tives of the characteristic function Xa) 


{(f^a ^ V®) X fl a }(J) a (x) — 2 Xa(*) 

ha 

(n a X V x )(/>a(x) = -2 --- Xa(x) = 


[a t • V a5 ]0 a (x) = -2a 


ha 
(at ‘ ha) 


-2 h~Xa{x ), 
aa 


t 7 

h a 


Xa{x) 


(C.67) 


(a • Vr)(6 • Vr) = (a n • Vji)(b n ■ Vr) + (a t • Vjt)(b n ■ Vr)+ 

+ («n • Vjt)(b t ■ Vr) + (a t • Vr)(6 ( • Vr) (C.68) 

(fln ' VR)(b n • Vr) = (tt n X Vr) • (f>n X Vr) + (<J n • b n )VR 

In the following we provide explicit expressions for matrix elements of the kernels G, P <3? 
and IT. 

Evaluation of matrix elements of G. 

f dS x dSy¥a(x)G(x,y)V b (y) = 

= ~ j dS x dSy(p a (x)(p b (y)[(a ■ b)gs(R) + jj(a ■ V R )(b ■ V R )gsc(R)] 

T a T b 5 

(C.69) 
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J dS x dS y cf> a (x)(f) b (y)(an • V H )(&t • V R )g S c(R) = ] dS x dS y M^ n • R)g'sc(R) 

T a T b 

J dS x dS y cf> a (x)(f) b (y)(a t • V R )(b n • V R )g S c(R) = 2( °^ K) J dS x dS y fo{y)(b n ■ R)tf sc (R) 

T a T b 

J dS x dS v Mx)Mv)(°t • Vr)(6* • V R )g SC (R ) = 4( °*' ^ J dS x dS y gsc(R ) 

T a T b 

j dSxdSy^a^x^cfrb^y^^cin • VR)(b n * V r)9sc(R') = 

T a T b 

= a n b n J dS x dSy{(l) a (x)(f) b (y)(n a -h b )[k 2 s gs{R)-kcgc{R)}-^-j^- i gsc{R)} 

TaT b 

(C.70) 

By collecting the above terms we obtain 


J dS x dS y (/)a(x)gsc(R)(a * Vr)(6 • V R )(j) b (y) = 

T a T h 

= a n b n {fi a • n 6 ) / - fec3C , (^)] + 

T a T h 

+ 2o J d S x dSy<l>a(x)(R-h a )^s C {R)^ 

T a T b 

H- ^——— J dS x dS y <fi b (y)(R • n b )g' sc (R)+ 

T a T b 

+ 4 (a-h a )(b-h b )-a n b n (l a -l b ) r dSxdSygsc{R) 
hah b J 

T a T b 


In deriving the above we used the following auxiliary relations 

n a xh a = l a 


For the linear basis function associated with node we have 

V x tl> a (*) = -^T 

J dSgsc{R)^ x 4 >a(x)] = -h a J dSgsc(R) 

For the constant basis function associated with triangle T’ and vertex v b we have 

3 

v a! ^( x) = - h b 5 b (r) 

6=1 
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[ dSgsc (R) [V x (pa (&)} = -ha [ dSgsc(R) 
J Jdl a 


Complete matrix element of G for the linear basis functions 
/ iS„dS y * b (y)G(R)* a ( x ) = (2-9/, + '*M>> ■£)(** 

TaT b 

4[(b ■ h b )(a • h a ) - (b • h b )(a • re a )(f a • l b )] 2(fe 6 • b)(a • re a ) 2(h a • a)(b- n & ) 

ykgh a h b 3 /j,kghb 4 l^kgh a 5 ’ 

(C.72) 

where 


h = J dS x J dSy4> a (x)(j) b (y)gs(R ), 

Ta T b 

h = j dS x JdS y <j> a (x)<j> b (y)[klgs(R) - k 2 c g c (R )\, 

T a T 6 

h = J dS x JdS y gsc(R ), 

T a T 6 

I 4 = J dS x J dSy(f> a (x)(n a -R)g'sc(R), 

T a T b 

h = J dS x J dS y (j) b (y)(n b ■ R)g' sc (R). 

Ta T h 


(C.73) 


Evaluation of matrix elements of P. 

J dS x dSy[* a (x ) • f(x, y) ■ * b (y)\ = J dS x dS y Mx)My)[bfa\, (C.74) 


where 

(bfa) = 


X 


-( n b ■ b)(a ■ V R )gc(R) + {(a ■ b)(h b ■ Vr) + (h b ■ a)(b • V R )}g s (fi)+ 


A -\- 2 /i 

+ -ro( a ' Vr)(& • Vfl)(n(, • VR)[gs(R) ~ gc(R)}, 


% 


(C.75) 


The operator P needs to be sandwiched between the piecewise linear basis function asso¬ 
ciated with vertex a of a triangle a and a piecewise constant function associated with triangle 
b : 

We will use the following decompositions of the V r into its tangential and normal compo¬ 
nents 


^R) — Q'nipa ’ T dt(ia ' 


(C.76) 
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Subsequently we evaluate the matrix elements 

J dS x dS v cf) a (x)cf) b (y)(n b ■ b))(a ■ R)g' c (R) = (fi b ■ b)[a n Gi + a t (t a • h a )G 2 ], 

(h b ■ a) J dS x dSy4) a (x)(j) b (y)(b ■ V R )g s (R) = (n b ■ a ) + b t (t b ■ h b )G A ), (C.77) 

(a • 6) J dS x dS y 4> a (x)(j) b (y)(n b ■ V R )g s (R) = (a • b)G 5 , 

The product of the three projections of the V# appearing in the matrix element of P can 
be suitably written as 

(h b ■ V)(6 • V)(a • V) = a n b n (h a • V)(n 6 • V)(n 6 • V)+ 

+ atbt(n b ■ V)(t a • V)(i& • V)+ 

■i - a v,bt(R b ■ • V)+ 

+ atb n (n b ■ V)(t a • V)(n& • V) = ^ 

= a n b n (h b • V)[(n a • n fe )V 2 - (n a x V) • (n 6 x V)]+ 

+ atbt(n b ■ V)(t a • V)(i& • V)+ 

+ a n b t (i b • V)[(n a • rift) V 2 - (n a x V) • (n 6 x V)]+ 

+ a t b n (i a ■ V)[(n 6 • n fe )V 2 - (n b x V) • (h b x V)] 

When sandwiched between the basis functions in the matrix element integrals the above terms 
yield 


{a n b n (n a ■ n b )4> a (x)<t> b {y) + 

+ 


^nbnija ' ^&) “1“ 


h a h b 


(n b ■ V)\gs(R) - gc(R)}+ 


0 J nbt(R'a ' rrfe)(i& ■ hh) ^ ( 7 ?) -|- ' fl'b)\t'a ' ^a) ^ 


+ 


h b 

{fi a ' rr&)(^a ' ^&) 
h a h b 


h n 


[k 2 s g s (R)-k 2 c gc(R)} + 


[a n b t (i b ■ h b )<f> a (R)5(r b ) + a t 6 n (i a • /r a )</>ft(i?)5(r a )][gf S ( J R) - S'cK-R)]- 


(C.79) 


Now the integrals of the term appearing in the above formula 

J dS x dS y 4> a {x)(f) b {y){fi b ■ V)(o • V)(& • V)bs(iJ) - <7c(#)] = 

7 / - - , /j —Q j nbn(j'a ’ P) T ’ h a )(tb ' ^ 

= a n b n (n a • n b ) G 6 + 4-—-G 7 + 

ha hfr 

CLnbti'fl'a ’ W'b)(ib ’ hb) ^ * h b )(t a • /l a ) 

+ ^;-Cj- 8 H-;-Cxg + 


(C.80) 


+ 4 


hb 

(h a • h b )(i a • Z 5 ) 
ha h b 


hn 


Q j nbt(i , b ■ T ^th n (£ a • ha)Gii]. 
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Complete matrix element of r for the linear basis functions 


J dS x dSy'P b (y)f(R)'J' a (x) = 

T a T b 

= ^ 2^ {’ b)[a n Gi + at(t a ■ h a )G 2 ] + (fi b ■ a)[b n Gs + bf(t b ■ h b )G^\ + (a • b)Gs}+ 

. 2fl n 6 n (fi a • fib) „ 8[ &nbn(Ja ' lb) “1“ ' ^a)(^6 ' ^fe)] . 

H-To-(j-0 H- - - - - —K -(-J7 + 


k 2 s 


+ 


+ 


h a h b k 2 s 

4a n b t {h a ■ n b )(t b ■ h b ) 4 a t b n {fi a ■ h b )(i a ■ h a ) 


h b kg 

&(ffa ' f^b) (Jo, ' lb) 


-G 8 + 


h a k 2 


g 9 + 


P’nbtO'b ' h b )G 10 “I - Clfb n (i a ■ fla)Gi±\. 


h a h b kg 

(C.81) 

The evaluation of the above matrix elements reduces to the calculation of the following 
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basic integrals 


Gi = j dS x j dSy4>a{x)4> b {y)(fi a ■ R)g' c (R), 
T a T b 

G 2 = J dS x J dS y 4> b (y)gc{R ), 

Ta T b 

Gz = J dS x J dSy{h a ■ R)g' s (R) = 

Ta T b 

G 4 = J ds x j dSygs(R), 


T a T b 


— 


Ga = 


G7 ~ 


J dS x J dSy<f> a {x)<f> h (y){fi h ■ R)g' s {R), 

Ta T b 

I ds x I dS v Ux)My)( n„ ■ mhs(R) - khc(R)], 

(C.82) 

Ta T b 

= J dS x J dSy(n b ■ R)[k 2 s gs(R ) - k 2 c g c (R)}, 

T a T b 

Gs = J dS x f dS y Mx)[k 2 3 gs(R ) - k 2 c g c (R)}, 

T a T b 

Gg = j dS x J dl y (/) b (x)[k 2 s gs{R) - kcgc(R)}, 

Ta T b 

Gio = J dS x J dl y (/)a(x)[gs(R) - g c (R)}, 

Ta l b 

Gn = hi dS y (f) b (x)[gs(R ) - gc(R)], 

la T b 

Evaluation of matrix elements of <P. Matrix elements of ^ Complete matrix element of 
for the linear basis functions 


J dS x dS y [V a (x ) • S(x,y) ■ * b {y )] = J 


= / dS x dS y 4> a {x)(t) b (y)[b®a ], (C.83) 


with 


A 


(b&a) = — - (h a ■ a)(b ■ V R )g c (R ) + {(a • b)(h a ■ Vr) + (n a ■ b)(a ■ 'V R )}gs(R)+ 


A H - *2/1 

+ Jr(a • V fl )(6 • V R )(n a • V r )[ 9s (R) - g c (R)], 
K s 


(C.84) 
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can obtained by transposition i.e., the following replacement in the expression for matrix 
elements of r 

a —y b 


/ 

TaT b 


dS x dS y \P b (y)<l>(x, yR)\P a (x) = 


A 


; {(n a • a)[b n Fi + b t (t b ■ h b )F 2 ] + (n a • b)[a n F 3 + a t (t b ■ h b )F 4 \ + (a • b)F 5 } + 


+ 


+ 


+ 


A -\- 2 /I 

2dnb n {'fla ’ 771 8 [ ttnbnija * lb) “ 1 “ ’ ^b)] 


kg 


-Fq + 


h a h b k 2 s 

4 a n b t (h a ■ h b )(t b ■ h b ) 4a t b n {fi a ■ fi b )(t a ■ h a ) 


Fj-\- 


h b k 2 s 

8(^a ’ 'f^b)O j a ' lb) 
h a h b kl 


-F$ + 


h a k 2 g 


Fq + 


P>nbt(ib ' b'b )^10 4 “ ' ^o)-^ll]- 


(C.85) 
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Fi= J dS x j dSy4> a {x)(j) b (y)(h b ■ R)g' c (R), 

T a T b 

F 2 = J dS x J dSy<j> a (x)g c (R), 

Ta T b 


Fa = 


F 3 = J dS x J dS y (n b • R)g's(R) = 

Ta T b 

F 4 = J dS x J dS y g s (R ), 

T a T b 

JdS x JdSy<j>a(x)(j>b(y)(na ■ R)g's(R), 

Ta T b 

= f dSx JdSyMx)My)(n a ■ R)[kgg s (R) - k 2 c g c (R )\, 

T a T b 

= J dS x J dSy(h a • R)[k 2 s gs(R) - k 2 c gc(R)}, 

T a T b 

J dS x J dS y (f> b (y)[k 2 s gs(R) - kcgc(R)}, 


Fq — 


F 7 = 


T a T b 

F 8 = 

Ta T b 

F 9 = j dS x J dl y (j> b (y)[kggs(R) - k 2 c g c (R )\, 

Ta T b 

F w = [ dS x [ dl y <j>a(x)\gs(R) - gc(R)], 


Ta lb 


F n = 


J dl x J dS y (f> b {y)\gs(R) - gc(R)], 


(C.86) 


la Tb 

Evaluation of matrix elements of W. Complete matrix element of W for the linear basis 
functions 


j dS x dS y [# a (x ) • W(x,y) • 9 b (y)] = J 


dS x dS y cf> a (x)My)[bWa}, (C.87) 
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with 


(b-W -a) = - A + ^ a n frnffcCft) + 

“t - ^ _|_ 2 ^ [®ra(^6 ’ V x){b ' V x) b n (fia ' V *) (® ‘ V aOl^C(-R) - ! - 

+ m{(« • 6)(«a • va.)(n 6 • Vj) + (n 6 • o)(n a • V x )(b • Va,)+ 

+ {h a ■ b)(a ■ V x )(n b ■ V x ) + (h a ■ n b )(a ■ V x )(b ■ V x )}g s (R)+ 

+ %(fi a • V*)(n 6 • V x )(a ■ V x )(b • V*)^) - g c (R)]- 
k s 

We note that 

Q'nifi'b ’ ^£c)(^ ' Vjc) + b n (fl a • V £C )(tt • Va;) = 

(^6 ^ ^ x) (^o ^ ^ x) ]H“ 

+ a n bt(fl b • V x ){ib ^ x) 4“ a tbn{ia ’ ^ x){fb a ’ ^x\ 


Subsequently we evaluate the integrals 

J dS x dSy<fi a ( K x)(fr b (y)[a n (Ti b • V x )(b * + b n (ft a • V £C )(tt • V £C )]^ , (y(i?) = 


T a T 6 


^nbn^jrT A&nbn 2 CL n bf{ib * ^6 T t r 2(2 tbnOa ’ h a 


- 2 k z c a n b n W 2 - 4 ^W 3 - 4 ^W 4 - —’~ a W 6 . 

^ h b K 


(C.88) 


(C.89) 


(C.90) 


|(a • 6 )(n a • V x )(Ti b ' 4 ~ (^6 * * ^x)(b * ^£c) 4 “ 

4“ (^a ' ^x){fbb ' ^£c) 4“ (^a ' ’ ^x)0 ’ ^£c)} = 

{(a • 6 )(n a • V^)(n a • V*) + b n (h b • a)(n a • V*)^ • V ffi ) + b t (n b • a)(n a • • V*)+ 

4“ CL n (fl a • b)(fl a • V x){fbb ^x) 4“ Q'tiP'a ' bk){t a ‘ V x){^b ’ ^x)~\~ 

^nbnij^a ‘ f^b)ix^a ‘ ^x){^b ^tc) 4“ ^nbtif^a * ^b){^a ’ ^cc)(^6 ‘ ^a?)4“ 

4“ CLtbnix^a ' Tib) ft a * ^ x){f^b ^ x) 4“ a tbt{^a ‘ ^ b)0a ‘ ^x)0b ‘ ^cc)} = 

"[((X • 6 ) 4“ ^ 71(^6 ' ^0 4h CL n {ri b • ft) 4~ Q j nbn(.'H / a ' * ^£c)(^fr * ^£c)4“ 

4“ 4“ Q'nbtifi'a ’ ^&)}(^a * ^ x)0b ’ ^cc)4“ 

+ {a t (n a • 6 ) + a t b n (h a • rx 6 )}(t a • V x ){ft b • V a3 )+ 

4“ {^^^(tXq, • tX5)}(t a • ^x)0b ’ V,). 

(C.91) 
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Finally 


{(d • b)(fl a • ^ x){fl b • Vgr;) + (n b * Q^)(fl a • V X )(b ' Va;) + 

"F (^a ’ ^x){Hb * Vas) “1“ (^a ' 'b'b)( a ' ^x)(b • V^)} = 

= {[(a • b) + b n (h b ■ a) + a n (n b ■ a) + a n b n (h a ■ n b )](n a ■ h b )}V 2 x + 

{[(a • 6) + b n (n b ■ a) + a n (n b ■ a) + a n b n (n a ■ h b )](h a ■ n b )}(n a x V*) • (n fe x V a! )+ (C.92) 
“1“ \ptiRb ‘ ®) “I - Q'nPtifi'a, ‘ ^b )}' ^x)(Pb ' 

+ {ot(^a • 6) + atb n (n a ■ Ti b )}(t a ■ ^x^iRb ■ ^x)4~ 

+ {o, b b b (n a ■ 'b.ft)}(ta * ^x)(ib ' v*). 

/ • b)(n a • V x )(n b ■ V x ) + (h b ■ a)(n a • V cc )(b • Va,)+ 

T a T b 

+ (n a • 6)(a • V x )(n b • V*) + (n a • n 6 )(a • Va,)(b • V*)} 5 s(i?) = 



W 8 ]+ 


+ 



2 [bt (^6 ' fl) “I - Q"nPtipb a • • fl b j ‘ 2 \o J i(fl a • b) + Q b b n (fi a • Tlb)](ia ' ^a) 

-;-”9-:- 


WlO+ 


+ 


(C.93) 


(»a • v a! )(n 6 • V x ){a ■ V x )(b ■ V x )\g s (R) - gc{R )] = 
Pa(x)cf) b (y)(n a ■ h b )(a • V x )(b ■ V x )[k 2 s g s (R) - k 2 c g c {R)] + 


(C.94) 



J dS x dS y cf> a (x)cf> b (y)(a • V R )(b • V R )feCR) - k 2 c g c (R)} 



(C.95) 


4[(a * fla)(b * /l^) C^n^nija ' ^&)] TJ r 

H--—:-n / i4 


J dS x dS y (a • V fl )(b • V R )fe(iJ) - g c (R)} = 


= a n b n (h a ■ n b )Wu + 2 a n (b t ■ h b )Wi 5 + 2 b n (a t ■ h a )Wie+ 
+ 4[(a • b. a )(b • fih) a n b n (l a • /ft)]VKi 7 . 
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Final expression for matrix element of W 


J dS x dS y Mx)My)(b ■ Wa) = 

TaT b 


Xkl Txr 

c W 2 + 


A H- 2 fl 

2 X/l r 2 7 T TT A Tjr A&vJOn JJT 2 CL n bf(tb ' ^6 yrr 2 Cifb n {t a ’ TT7 'I . 

+ T . ^ {-2/^a n ft n VF 2 - 4-^-VF 3 - 4—^ 2 ~W 4 -^-VF 5 ,-£-VF 6 }+ 


A -f- 2 fi 


h 2 

a b 


h b 


h n 


+ /i{[(a • b) + b n (n b ■ a) + a n (n b ■ a) + a n b n (h a ■ h b )](n a ■ n b )[-k 2 s W 7 - ^—^W 8 }+ 
-2 [b t (n b ■ a) + a n b t (h a ■ h b )](t b ■ h b ) 2 [a t (n a ■ b ) + a t b n (n a ■ n b )]{i a ■ K) 

H-;-W 9---tVio+ 


h b 

_l_ ^tbt{fl a * flb){t a * h a ){tb • hb) 


h n 


hah b 


4/i(n a • h b ) 


2 CL n (ht * h'b) TT7 - 1 2b n (ci t • h a ) 


c a ,v b) r 7 / - - \ T i/ , , 

T2-\ a n^n( n a * n 6jW / ll H-7-^12 + 

tin 


hb 


W 13 + 


_l_ 4[(fl • ha){b • hb) a n b n (l a • lb)] 


4(ia * lb) 
ha h b 


h a h b 

{a n b n (h a • h b )W u + 


+ 2 a n {pt • fo&)VFi 5 + 2 b n [at • /&a)^i6 + 4[(a • h a )(h • /&&) — a n b n (l a • /&)]VFi 7 }. 


(C.97) 
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where 


W 2 = J dS x dSy<t> a (x)<t> h {y)gc{R ), 

TaT b 

w 3 = f dl a J dSyMv)9c(R), 

la Tb 

Wi = J dS x J dl b (f> a (x)g c {R)+ 

T a lb 

w 5 = J dS x dSy(f) a (x)(h b ■ V x )gc(R). 

T a T b 

W e = dS x dS v (f> b (y)(n a ■ V x )g c (R), 

T a T b 

W 7 = J dS x dSy(f> a (x)(f> b (y)gs{R), 

T a T b 

W 8 = f dS x dS y gs(R), 

T a T b 

Wg = j dSxdSyfia^x^fib • V x )gs(R)j 

TaTb 

= j dS x dSy(f) b (y)(n a • V a3 )^ f 5 '(-R), 

T a T h 

Wu= J dS x dSyMx)<t>b(y)[ksgs(R) - k 4 c gc(R)}, 

TaTb 

w 12 = J dS x dSyMx)(R-n a )[k 2 s g' s (R)-k 2 c g' c (R)}, 

T a T h 

W 13 = j dS x dS y (f> b (y){R ■ h b )[k 2 s g' s {R) - k 2 c g' c (R)\, 
T a T h 

W 14 = J dS x dS y [k 2 s gs(R) - k 2 c g c (R)\, 

T a T h 

^ 15 = J j dS x dl b {R-h a )\g' s {R)-g' c (R)l 

T a lb 

W 16 = j j dl a dS y (R • nMiR) - k 2 c g' c (R)}, 

la Tb 

W 17 = J dl a dl b [k 2 s g s (R) - k 2 c g c (R)}. 


(C.98) 
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C.5.4 Behavior of the kernels near the Green function singularity R — 0 

We note that the displacement Green functions G(R) and and the surface stress tensor Green 
function r(r,r') appearing in the surface integral equation exhibit an important property in 
the limit R —> 0. Since the function 


gs(R) ~ gc(R) 


\ /gi k s R _ I e i k c R_l\ 

4jt V R R ) 


(C.99) 


is regular for R —>> 0 (due to cancellation of the singularities in the two Green functions), the 
second terms of the Green functions G(R) and n • U(r,r') are also nonsingular, while, without 
the cancellation, they would have contained a 1/i? 3 singularity. The reduced degree of singu¬ 
larity is particularly important in the discretization of surface integral equations (Sec. C.4). 
Specifically we have for small R 

[ gs(R ) - gc(R)] ->• -7- {i(^5 - kc) + 5[(Vs) 2 - (i kc) 2 ]R + •••} 


as well as 


(n ■ V„)te(i?) - sc m = (n • R)±-±\g s (R) - gc (K)] -+—(»• A)[(*l - k 2 c ) + ...] 

are finite at R = 0 This fact facilitates discretization of the integral equations and computation 
of the matrix elements. 
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